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

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 datos con errores distribuidos normalmente. Este es el modelo de regresión más 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 se relaja, introduciendo la función de enlace. En segundo lugar, las distribuciones de errores distintas de las normales se pueden modelar

Modelos Lineales Generalizados

Un modelo de regresión define la distribución de una variable de respuesta (a menudo denotada 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 utilizado, la regresión lineal ordinaria, modela 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 se puede representar 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 monotónica 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 "link". Los ejemplos incluyen el link logit (sigmoid) y el link del registro. Además, y puede tener una distribución no normal, como el binomio o Poisson. Por ejemplo, una regresión de Poisson con vínculo de registro y un único predictor 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 coches de varios pesos que no superan una prueba de kilometraje. Los datos incluyen observaciones de peso, número de coches probados y número fallado.

% 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 gráfica de la proporción de coches que fallan, en función del peso. Es razonable suponer que los recuentos de errores provienen de una distribución binomial, con un parámetro de probabilidad P que aumenta con el peso. Pero, ¿cómo debe depender exactamente P del 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 limitadas. Esto infringe una de las suposiciones necesarias para ajustar un modelo de regresión lineal simple.

El uso de un polinomio de orden superior puede parecer para 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 es superior a 4000; de hecho, se volverá negativo para valores de peso más grandes. Y, por supuesto, la suposición de una distribución normal sigue siendo violada.

En su lugar, un mejor enfoque es usar para ajustar 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 pronosticadas al rango [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 asimptotas a cero y una como peso se vuelve pequeña o grande.

Diagnóstico de modelos

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 de dos modelos para determinar si un término cuadrado mejoraría significativamente el ajuste.

[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 gráfica 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 residuos 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 diferentes varianzas.)

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

Evaluación de las predicciones del modelo

Una vez que estamos satisfechos con el modelo, podemos usarlo para hacer predicciones, incluyendo límites de confianza de computación. Aquí predecimos el número esperado de coches, de 100 probados, que fallaría 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 (predeterminada).glmfit Para la distribución binomial, el enlace canónico es el logit. Sin embargo, también hay otros tres enlaces que son sensibles 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 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.