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.

Datos de ajuste con modelos lineales generalizados

En este ejemplo se muestra cómo ajustar y evaluar modelos lineales generalizados mediante y.glmfitglmval La regresión lineal ordinaria se puede utilizar para ajustar una línea recta, o cualquier función que sea lineal en sus parámetros, a los datos con errores normalmente distribuidos. Este es el modelo de regresión más comúnmente utilizado; sin embargo, no siempre es realista. Los modelos lineales generalizados amplían el modelo lineal de dos maneras. En primer lugar, la asunción de linealidad en los parámetros es relajada, introduciendo la función de enlace. En segundo lugar, las distribuciones de errores que no sean normales se pueden modelar

Modelos lineales generalizados

Un modelo de regresión define la distribución de una variable de respuesta (a menudo denotado genéricamente como y) en términos de una o más variables predictoras (a menudo denotadas x1, x2, etc.). El modelo de regresión más comúnmente utilizado, la regresión lineal ordinaria, modelos y como una variable aleatoria normal, cuya media es la función lineal de los predictores, B0 + B1 * x1 +..., y cuya varianza es constante. En el caso más simple de un solo predictor x, el modelo puede representarse como una línea recta con distribuciones Gaussianas sobre cada punto.

mu = @(x) -1.9+.23*x; x = 5:.1:15; yhat = mu(x); dy = -3.5:.1:3.5; sz = size(dy); k = (length(dy)+1)/2; x1 =  7*ones(sz); y1 = mu(x1)+dy; z1 = normpdf(y1,mu(x1),1); x2 = 10*ones(sz); y2 = mu(x2)+dy; z2 = normpdf(y2,mu(x2),1); x3 = 13*ones(sz); y3 = mu(x3)+dy; z3 = normpdf(y3,mu(x3),1); plot3(x,yhat,zeros(size(x)),'b-', ...       x1,y1,z1,'r-', x1([k k]),y1([k k]),[0 z1(k)],'r:', ...       x2,y2,z2,'r-', x2([k k]),y2([k k]),[0 z2(k)],'r:', ...       x3,y3,z3,'r-', x3([k k]),y3([k k]),[0 z3(k)],'r:'); zlim([0 1]); xlabel('X'); ylabel('Y'); zlabel('Probability density'); grid on; view([-45 45]); 

En un modelo lineal generalizado, la media de la respuesta se modela como una transformación no lineal monótona de una función lineal de los predictores, g (b0 + B1 * x1 +...). La inversa de la transformación g se conoce como la función de "enlace". Algunos ejemplos son el vínculo logit (Sigmoid) y el vínculo de registro. Además, y puede tener una distribución no normal, como el binomio o Poisson. Por ejemplo, una regresión de Poisson con enlace de registro y un predictor único x se pueden representar como una curva exponencial con distribuciones de Poisson sobre cada punto.

mu = @(x) exp(-1.9+.23*x); x = 5:.1:15; yhat = mu(x); x1 =  7*ones(1,5);  y1 = 0:4; z1 = poisspdf(y1,mu(x1)); x2 = 10*ones(1,7); y2 = 0:6; z2 = poisspdf(y2,mu(x2)); x3 = 13*ones(1,9); y3 = 0:8; z3 = poisspdf(y3,mu(x3)); plot3(x,yhat,zeros(size(x)),'b-', ...       [x1; x1],[y1; y1],[z1; zeros(size(y1))],'r-', x1,y1,z1,'r.', ...       [x2; x2],[y2; y2],[z2; zeros(size(y2))],'r-', x2,y2,z2,'r.', ...       [x3; x3],[y3; y3],[z3; zeros(size(y3))],'r-', x3,y3,z3,'r.'); zlim([0 1]); xlabel('X'); ylabel('Y'); zlabel('Probability'); grid on; view([-45 45]); 

Ajuste de una regresión logística

Este ejemplo implica un experimento para ayudar a modelar la proporción de automóviles de varios pesos que fallan una prueba de kilometraje. Los datos incluyen observaciones de peso, número de coches probados y número fallido.

% A set of car weights weight = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]'; % The number of cars tested at each weight tested = [48 42 31 34 31 21 23 23 21 16 17 21]'; % The number of cars failing the test at each weight failed = [1 2 0 3 8 8 14 17 19 15 17 21]'; % The proportion of cars failing for each weight proportion = failed ./ tested;  plot(weight,proportion,'s') xlabel('Weight'); ylabel('Proportion'); 

Este gráfico es una trama de la proporción de coches que fallan, en función del peso. Es razonable suponer que los recuentos de fallas provienen de una distribución binomial, con un parámetro de probabilidad P que aumenta con el peso. Pero ¿cómo exactamente debe depender P de peso?

Podemos intentar ajustar una línea recta a estos datos.

linearCoef = polyfit(weight,proportion,1); linearFit = polyval(linearCoef,weight); plot(weight,proportion,'s', weight,linearFit,'r-', [2000 4500],[0 0],'k:', [2000 4500],[1 1],'k:') xlabel('Weight'); ylabel('Proportion'); 

Hay dos problemas con este ajuste lineal:

1) la línea predice proporciones menores que 0 y mayores que 1.

2) las proporciones no se distribuyen normalmente, ya que están necesariamente delimitadas. Esto infringe uno de los supuestos necesarios para ajustar un modelo de regresión lineal simple.

El uso de un polinomio de orden superior puede parecer ayudar.

[cubicCoef,stats,ctr] = polyfit(weight,proportion,3); cubicFit = polyval(cubicCoef,weight,[],ctr); plot(weight,proportion,'s', weight,cubicFit,'r-', [2000 4500],[0 0],'k:', [2000 4500],[1 1],'k:') xlabel('Weight'); ylabel('Proportion'); 

Sin embargo, este ajuste todavía tiene problemas similares. El gráfico muestra que la proporción ajustada comienza a disminuir a medida que el peso va por encima de 4000; de hecho, se convertirá en negativo para los valores de peso más grandes. Y, por supuesto, se sigue violando la suposición de una distribución normal.

En su lugar, un mejor enfoque es usar para ajustarse a un modelo de regresión logística.glmfit La regresión logística es un caso especial de un modelo lineal generalizado y es más adecuada que una regresión lineal para estos datos, por dos razones. En primer lugar, utiliza un método de ajuste que es adecuado para la distribución binomial. En segundo lugar, el enlace logístico limita las proporciones previstas al intervalo [0,1].

Para la regresión logística, especificamos la matriz predictora y una matriz con una columna que contiene los recuentos de errores y una columna que contiene el número probado. También especificamos la distribución binomial y el enlace logit.

[logitCoef,dev] = glmfit(weight,[failed tested],'binomial','logit'); logitFit = glmval(logitCoef,weight,'logit'); plot(weight,proportion,'bs', weight,logitFit,'r-'); xlabel('Weight'); ylabel('Proportion'); 

Como indica esta gráfica, las proporciones ajustadas asíntota a cero y una como peso se vuelve pequeña o grande.

Diagnósticos de modelo

La función proporciona una serie de salidas para examinar el ajuste y probar el modelo.glmfit Por ejemplo, podemos comparar los valores de desviación para dos modelos para determinar si un término cuadrado mejoraría el ajuste significativamente.

[logitCoef2,dev2] = glmfit([weight weight.^2],[failed tested],'binomial','logit'); pval = 1 - chi2cdf(dev-dev2,1) 
 pval =      0.4019  

El valor p grande indica que, para estos datos, un término cuadrático no mejora significativamente el ajuste. Una trama de los dos ajustes muestra que hay poca diferencia en los ajustes.

logitFit2 = glmval(logitCoef2,[weight weight.^2],'logit'); plot(weight,proportion,'bs', weight,logitFit,'r-', weight,logitFit2,'g-'); legend('Data','Linear Terms','Linear and Quadratic Terms','Location','northwest'); 

Para comprobar la bondad del ajuste, también podemos ver una gráfica de probabilidad de los residuales de Pearson. Estos se normalizan de modo que cuando el modelo es un ajuste razonable a los datos, tienen aproximadamente una distribución normal estándar. (Sin esta estandarización, los residuos tendrían varianzas diferentes.)

[logitCoef,dev,stats] = glmfit(weight,[failed tested],'binomial','logit'); normplot(stats.residp); 

La gráfica residual muestra un buen acuerdo con la distribución normal.

Evaluar las predicciones del modelo

Una vez que estemos satisfechos con el modelo, podemos utilizarlo para hacer predicciones, incluyendo la computación de los límites de confianza. Aquí predecimos el número esperado de coches, de 100 probados, que fallarían la prueba de kilometraje en cada uno de los cuatro pesos.

weightPred = 2500:500:4000; [failedPred,dlo,dhi] = glmval(logitCoef,weightPred,'logit',stats,.95,100); errorbar(weightPred,failedPred,dlo,dhi,':'); 

Funciones de enlace para modelos binomiales

Para cada una de las cinco distribuciones que admite, hay una función de vínculo canónico (predeterminado).glmfit Para la distribución binomial, el enlace canónico es el logit. Sin embargo, también hay otros tres enlaces que son sensatos para los modelos binomiales. Los cuatro mantienen la respuesta media en el intervalo [0,1].

eta = -5:.1:5; plot(eta,1 ./ (1 + exp(-eta)),'-', eta,normcdf(eta), '-', ...      eta,1 - exp(-exp(eta)),'-', eta,exp(-exp(eta)),'-'); xlabel('Linear function of predictors'); ylabel('Predicted mean response'); legend('logit','probit','complementary log-log','log-log','location','east'); 

Por ejemplo, podemos comparar un ajuste con el enlace de probit a uno con el enlace logit.

probitCoef = glmfit(weight,[failed tested],'binomial','probit'); probitFit = glmval(probitCoef,weight,'probit'); plot(weight,proportion,'bs', weight,logitFit,'r-', weight,probitFit,'g-'); legend('Data','Logit model','Probit model','Location','northwest'); 

A menudo es difícil para los datos distinguir entre estas cuatro funciones de enlace, y a menudo se hace una elección por motivos teóricos.