El uso de modelos lineales para predecir los resultados
Los ejercicios de laboratorio se centrarán en la regresión lineal simple y presentar e introducir la verosimilud como la base para evaluar el apoyo para los modelos alternativos (hipótesis). Los estudiantes también van a prender cómo usar la suma de cuadrados para estimar la incertidumbre de las predicciones del modelo. Hay enfoque en la interpretación (comprensión) de los valores estimados de los parámetros, la evaluación de las pruebas entre los modelos alternativos, y la predicción. El laboratorio también requerirá estudiantes ajustar los datos a 2-3 modelos de regresión lineal, comparar sus niveles de apoyo, interpretar los valores estimados de los parámetros, y utilizar los estimadors para predecir respuestas.
Las ideas generales de un modelo lineal
Un modelo lineal es simplemente un modelo en el que una respuesta (escrito frecuentemente como Y) se modela como una función lineal de uno o más predictores (factores explicativos, escrito frecuentemente como X) y los parámetros o coeficientes, escrito como b o, a veces la letra griega β. Los modelos lineales incluyen un amplio clase de modelos estadísticos incluido el análisis de la varianza (ANOVA), análisis de covarianza (ANCOVA) y modelos lineales generalizados (GLIM). Vamos a enfocar en los modelos de regresión en que los predictores son típicamente medidas continuas en lugar de niveles de los factores (como en ANOVA).
En el modelo de regresión lineal simple hay un solo predictor X y ambos el predictor y la respuesta Y se asume generalmente para ser en escala continua (real), aunque esto puede implicar transformaciones de datos que son originalmente en el recuento, binario, o de otro tipo no continuo, los bienes
escalas. La forma matemática básica de un modelo de regresión lineal simple es
Y i ^ = b 0 + b 1 X i
donde el "^" en la Y indica que la observación Y se está predicha por el modelo. En este modelo, hay 2 parámetros: la intersección b0 y la pendiente b1 . Como veremos más adelante, también hay un tercero parámetro representando lo bien que el modelo "se ajusta" a los datos. La interpretación de los parámetros de regresión simple es sencillo: el parámetro b0 representa la intersección, que es el valor que la respuesta Y toma cuando X = 0. La parámetro b1, la pendiente, representa el cambio en Y para cada unidad de cambio en X. Por el contrario, si tenemos una descripción verbal de estas relaciones es fácil convertirlo a una forma matemática. Por ejemplo, si yo digo que la respuesta es 3 cuando X = 0 y se doble cada unidad en aumento de X, puedo convertir esta afirmación en un modelo lineal como
Y i ^ = 3 + 2 X i
Ejemplo - Relación entre talla y peso en los peces del mar.
Nuestro ejemplo motivador es la relación entre la longitud en cm y el "peso" (masa en kg) de pescado de mar capturados en las redes de arrastre de fondo (NOAA 2003). Para crear los modelos lineales, los investigadores primera transformaron los datos mediante logaritmos naturales, de manera que el modelo real fue
log (M) = b0 b1 + log (L)
La interpretación de los parámetros de la escala de registro es que la intersección representa el valor de log (masa) cuando log (longitud) = 0, o cuando L = 1 cm; la pendiente es la tasa de aumento en el registro escalar por unidad de aumento de log (longitud). Los investigadores desarrollan modelos para varias especies de peces. Por ejemplo, el modelo de cazón suave masculina (Mustelus canis) basado en> 800 observaciones (pares de longitud de masas), en un rango de 40 a 135 cm fue
log (M) = 12.3265 + 2.9436log (L)
El guión R ilustra estos ejemplos y realiza gráficas.
La estimación de los parámetros de los modelos lineales
El problema básico de la estimación en la regresión lineal (y modelos lineales en general) es similar al el problema general de estimación que cubre en el último laboratorio. Es decir, tenemos un modelo estadístico donde la hipótesis se describen los datos, y ahora tenemos en la mano los datos (en este caso, pares x, y). Nuestras tareas es utilizar el modelo y los datos para derivar valores (estimaciones) de los parámetros que son en algún sentido óptima; estamos a punto de definir lo que entendemos por "óptimo".
Volvamos a nuestro modelo de regresión lineal simple en forma genérica. Estamos utilizando un modelo de la forma
Y = b0 + B1X
para explicar la variación en la respuesta i Y (en la regresión que convencionalmente suponemos que los valores de i X no varían al azar, sino que están en lugar fijo y mide exactamente; obviamente, esto es, en el mejor
una aproximación en la realidad). La diferencia entre la respuesta real y la predicha la respuesta se conoce como el residuo
ε i = Y i - Y i ^
o equivalentemente
Y = b0 + b1 X + ε i
Claramente, ε i representa cómo de cerca el valor predicho de Y para un valor dado de X en los partidos de la valor real de Y para ese mismo valor X. Si todos los ε i son iguales a cero, entonces el ajuste es perfecto (las predicciones coincidir exactamente con los datos). Nuestro enfoque de estimación será en cierta medida, dependerá de lo que asumir sobre la distribución estadística de las ε i, que a continuación, determinar el método más apropiado para la estimación de parámetros para optimizar el ajuste del modelo a los datos. Los 2 métodos básicos vamos a tener en cuenta son la estimación de mínimos cuadrados y la estimación de máxima verosimilitud.
Vamos a utilizar una muestra de 872 mediciones de talla-peso cazón para motivar a cada enfoque. Los datos se encuentran aquí. En primer lugar vamos a leer los datos y representar los datos en el original y registro
escalas transformadas-
> #read in length mass data
> dogfish<-read.csv("length-mass.csv")
> with(dogfish,plot(L,M))
> with(dogfish,plot(log(L),log(M)))
Sin duda, una relación lineal en la escala de registro parece bastante razonable para este tipo de datos. Nos pondremos en contacto un poco más específico acerca de lo que entendemos por "razonable" a medida que avancemos.
Estimación de mínimos cuadrados
Estimación de mínimos cuadrados hace lo que su nombre indica-que encuentra valores de los parámetros que minimizan la suma de los cuadrados de los residuos. Es decir, se trata de encontrar los parámetros estimados que minimizan la suma de los residuos al cuadrado
Σ (Y i - Y i ^) 2
La función R que peform estimación de mínimos cuadrados es la función () lm. Podemos ilustrar con este ejemplo de datos.
Los cálculos idénticos se pueden realizar en el LM R () del objeto.
El código
> #read in length mass data
> dogfish<-read.csv("length-mass.csv")
> with(dogfish,plot(L,M))
> with(dogfish,plot(log(L),log(M)))
> lm function in R
> reg<-lm(log(M)~log(L),data=dogfish)
> summary(reg)
calcula los estimados de mínimos cuadrados y superpone una gráfica de la ecuación de predicción resultante en los datos. También se produce una tabla de resumen. Aquí usted puede notar que los valores de p, así como
se producen las estadísticas de probabilidad. Pero, en realidad, no hemos hecho ninguna suposición sobre el distribución de los datos (no es necesario para los mínimos cuadrados), por lo que ¿de dónde salieron? Resulta el de los modelos lineales, los parámetros estimados por mínimos cuadrados y máxima verosimilitud estimaciones son idénticos, si se supone que los residuos se distribuyen normalmente. Es decir estamos
suponiendo que
ε i ~ N (0, σ)
es decir que estamos suponiendo que el modelo es imparcial (que significa que los residuos deben promediar cero) y tienen una varianza de σ 2. Más adelante veremos exactamente por qué esta relación funciona para
Errores distribuidos normalmente (pero no necesariamente para otras distribuciones). Si queremos invocar supuestos normales (o examinar la validez del modelo normal) podemos utilice el comando plot (reg) para trazar el objeto de modelo lineal. Esto producirá una serie de gráficos de diagnóstico que se pueden utilizar para diagnosticar problemas con el modelo (como la no linealidad) o violaciónes de los supuestos normales. Los comandos
> #diagnostic plots
> par(mfrow=c(2,2))
> plot(reg)
> savePlot(filename="Rbuilt_in_plots_reg",type="png",device=dev.cur())
> par(mfrow=c(1,1))
>
crear una serie de 4 parcelas de diagnóstico y guardarlos en un archivo en una sola página. Para 3 de los 4 estamos en busca de una dispersión aleatoria de los residuos; para las 4 ª (el QQ normal) que estamos buscando
datos que caen cerca de la línea recta. Un rápido examen de los gráficos de residuos y cuantiles sugiere que el modelo y los datos se comportan bien; Sin embargo, las pruebas de diagnóstico más formales son
incorporado en R.
Bondad de ajuste del modelo
Los modelos de mínimos cuadrados estimaciones ofrecen, por definición, el mejor modelo lineal capaz de que describe los datos que hemos utilizado para construir el modelo, y la varianza residual estimado es nuestra estimación de lo que está "de sobra." Podemos usar la función () lm derivar alguna otra estadísticas que son útiles para resumir la regresión. En primer lugar, si volvemos a la suma residual de cuadrados (la cantidad que es, por definición, minimizada bajo mínimos cuadrados) que tenemos en este
ejemplo, un valor de 5.940393
R ofrece esta directamente del reg <-lm (log (M) ~ log (L)) objeto a través de la desviación que se proporciona por deviance ()
>rss<- deviance(reg)
[1] 5.940393
Podemos comparar este valor a un modelo en el que no tenemos predictores, en otras palabras, el
interceptar único modelo, para obtener una estimación de las sumas totales de cuadrados variación de Y de su propia
significar. R se pone este valor a través de instalación de un único modelo de intercepción
> #Suma total de cuadrados (model de solo interceccion)
> ss_total <-with(dogfish,deviance(lm (log (M) ~ 1)))
> ss_total
[1] 803.8997
Por sustracción de la desviación residual de la desviación total que obtenemos la suma de cuadrados
explicada por el modelo de regresión.
> # Suma de cuadrados de regresión
> ssr <-ss_total-rss
> ssr
[1] 797.9593
Por último, podemos obtener la proporción de la desviación total del explicada por el modelo
> R2 <-ssr/ss_total
> R2
[1] 0.9926105
R2 (también conocida como el coeficiente de determinación) es un valor entre 0 (sin variación explicada
por el modelo) a 1 (100% de la variación explicada por el modelo) al 1 (100% de la variación explicada por el modelo). Es común pensar en R 2 en términos de "más alto es mejor". Sin embargo, como veremos 2 R se puede hacer arbitrariamente alto añadiendo términos al modelo de regresión, que en algún momento no se apoya en los datos. Esto se convierte en una cuestión bastante seria como la complejidad del modelo se eleva con respecto al tamaño de la muestra. Pensar en lo que significa, por ejemplo, para ajustar una recta a 2 puntos de datos (pista: sólo hay una manera de hacerlo!).
Evaluación de la incertidumbre en los modelos lineales
Cuando hemos ajustado un modelo lineal, hemos utilizado los datos (presumiblemente de un experimento o una muestra aleatoria) para estimar los parámetros y para hacer predicciones. Debido a que el modelo se basa en una muestra, contendrá incertidumbre inherente. Por lo tanto, una práctica habitual es proporcionar valores estimados de la variación o la confianza en los resultados del modelo. Hay 2 tipos básicos de incertidumbre que nos queremos evaluar e informar, y que están relacionados. Uno de ellos es la incertidumbre en los valores de los parámetros a sí mismos, es decir, en nuestras valores estimados de b0 y b1. Bajo estimación de mínimos cuadrados podemos obtener valores estimados de los errores estándares de los parámetros, que en última instancia se refieren de nuevo a lo bien que el modelo se ajusta a los datos (medida por RSS). Cuando asumimos Normalidad u otra distribución también podemos construido intervalos de confianza y realizar pruebas estadísticas basadas en la probabilidad.
El otro tipo de incertidumbre que nos interesa --- posiblemente más importante-es la capacidad del modelo para predecir y, en particular, su capacidad para predecir las afueras de los datos que utiliza para estimar los parámetros. Vamos a introducir la predicción de abajo, y volver al tema de la predicción y validación posterior en el laboratorio.
Predicción de los modelos lineales
Por encima nos hemos centrado en la descripción de lo bien que el modelo funciona para explicar la variación en los puntos (Y, X) de datos utilizado para construir el modelo. Por lo general, vamos a estar interesados en utilizar el modelo para la predicción, donde se observó un valor de X y (en función del modelo) predecir lo que el correspondiente valor de Y debe ser. Hay 2 tipos diferentes de las predicciones.
En la primera, que se remonta a los mismos valores de X (digamos, en otra muestra) y la predicción de un nuevo valor de Y. En el segundo, estamos usando un valor de X no en los datos originales, y quizás incluso fuera del rango de X que hemos observado en nuestra muestra. Intuitivamente, debemos hacer mejor en la predicción en el primer caso que en el segundo, y de hecho la predicción será la más precisa cerca de la media de los valores de X utilizado para ajustar el modelo. Otra distinción es si el interés se centra en la predicción de la media respuesta Y para un valor dado de X, o en la predicción de un valor particular de y que aún no ha sido observado. Estas distinciones se vuelven importantes para el cálculo adecuado de intervalos de confianza (que por ahora se producen bajo suposiciones normales), como se ilustra por los ejemplos a continuación.
La predicción se puede realizar ya sea mediante la obtención de los valores estimados y conectarlos (junto con valores de X) en la ecuación de regresión.
> #PREDICTION
> #predict Y at the values of length of 40, 60,80, and 140 cm
> #direct method
> b0<-coef(reg)[1]
> b1<-coef(reg)[2]
> new <- data.frame(L = c(40, 60,80,140))
> y_pred<-b0+b1*log(new$L)
> y_pred
[1] -1.4621266 -0.2734437 0.5699403 2.2105395
> mass_pred<-exp(y_pred)
> mass_pred
[1] 0.2317429 0.7607552 1.7681614 9.1206353
> y_pred<-predict(reg,new)
> y_pred
1 2 3 4
-1.4621266 -0.2734437 0.5699403 2.2105395
> mass_pred<-exp(y_pred)
> mass_pred
1 2 3 4
0.2317429 0.7607552 1.7681614 9.1206353
Como alternativa, se puede utilizar la función de predict() en R para producir predicciones, errores típicos, y los intervalos de confianza. Bajo supuestos normales, también podemos producir las parcelas con límites de confianza. Hay 2 clases de límites de confianza, en función de si estamos prediciendo la respuesta media o un individuo valor futuro. Bajo supuestos normales, también podemos producir las parcelas con límites de confianza. Hay 2 clases de límites de confianza, en función de si estamos prediciendo la respuesta media o un individuo valor futuro. Código para cada momento, para el ejemplo de longitud masa es
> #USING PREDICT FUNCTION
> #create normal confidence and prediction intervals and plot them
> #predict mean response
> ypred_c<-predict(reg,interval="c")
> #predict an unobserved individual responseS
> ypred_p<-predict(reg,interval="p")
Warning message:
In predict.lm(reg, interval = "p") :
predictions on current data refer to _future_ responses
> with(dogfish,plot(log(L),log(M)))
> matlines(log(dogfish$L),ypred_p,col="red",lty=c(1,2,2))
> ypred<-predict(reg)
>
> # Crear intervalos de confianza y predicción normales y crear gráficas
> ypred_c <-predict(reg, interval = "c")
> with(dogfish,plot (log (L), log (M)))
> with(dogfish,matlines(log (L), ypred_c, col = "red", LTY =c(1,2,2)))
>
Por defecto, las predicciones y los intervalos se calculan y se representan en el rango de la
datos X observados. Si queremos especificar un conjunto diferente de datos X podemos hacerlo, por ejemplo:
> new <- data.frame(L = c(40, 60,80,140))
> ypred2<-predict(reg, new, interval="p")
> ypred2
fit lwr upr
1 -1.4621266 -1.6249016 -1.2993516
2 -0.2734437 -0.4358264 -0.1110610
3 0.5699403 0.4076624 0.7322182
4 2.2105395 2.0480517 2.3730272
La estimación de máxima verosimilitud y modelos lineales generalizados
Como hemos visto anteriormente, la estimación de máxima verosimilitud consiste en la especificación del modelo (una distribución estadística). Los datos se tratan como constante, y valores de los parámetros se seleccionan para maximizar la función de verosimilitud. La única diferencia ahora es que bajo modelos lineales se imponen la restricción es que en lugar de encontrar los valores de los parámetros para la distribución original, estamos encontrando valores correspondiente al modelo lineal. En concreto, nuestro modelo está ahora condicionado al valor de X a través del modelo lineal. El examen de esta probabilidad revela que la maximización es equivalente matemáticamente para minimizar los residuos del modelo. Por lo tanto, si suponemos una probabilidad normal de los mínimos cuadrados y las estimaciones de máxima verosimilitud son intercambiables.
Como se mencionó anteriormente, la función lm () proporcionó ambos los estimados de mínimos cuadrados y los MLEs (Y valores asociado de la log-verosimilitud) bajo supuestos normales. Un procedimiento más general para permitir la especificación de otras probabilidades se proporciona mediante modelos lineales generalizados (glm) y la función glm() en R. Los modelos lineales generalizados pueden ser apropiados si los supuestos normales ya no aplican (por ejemplo, con el recuento o respuestas binarias) o si el examen de los residuos u otras pruebas muestran desviación de la normalidad. El procedimiento glm() permite tanto especificación de una distribución de la "familia" (por ejemplo, normal o de Gauss, Poisson, binomial) y una transformación o función de enlace (identidad [por defecto], registro, logit, etc.).
Como ejemplo, vamos a considerar un estudio de nidos curruca a 150 m de 100 a 1500 m en el sur de las montañas Apalaches. En cada elevación, entre n = se encontraron 30 y 40 nidos y se los seguin hasta que el éxito (incipiente) o fracaso (destrucción de nidos, huevos o polluelos). La respuesta era por lo tanto el número de nidos y éxito en n ensayos en cada elevación. Los datos para el estudio son aquí y el código de R esta aquí. Podemos construir un modelo lineal (en la escala logit) con la asunción de error binomial para éstos datos como
> rm(list=ls())
> data.dir<-"C:/Users/Mike/Dropbox/company/fulbright/web/2"
> setwd(data.dir)
> nest_data <-read.csv ("nest_success2.csv")
> glm1 <-glm (y / n ~ elev, weight = n, family = binomial (link = logit), data = nest_data)
> summary (glm1)
Call:
glm(formula = y/n ~ elev, family = binomial(link = logit), data = nest_data,
weights = n)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4313 -0.4732 0.0273 0.8200 1.9239
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.6365905 0.5628019 -8.238 <2e-16 ***
elev 0.0087353 0.0009975 8.757 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 290.841 on 9 degrees of freedom
Residual deviance: 10.255 on 8 degrees of freedom
AIC: 30.857
Number of Fisher Scoring iterations: 5
La generación de los intervalos de confianza para predicciones con glm () es un poco más complicado de lo debido al hecho de que no estamos asumiendo una única distribución. Aquí vamos a utilizar una aproximación normal
de las predicciones en la escala logit (que deberían por aproximadamente cierto dados máxima verosimilitud) junto con la predicción y su error estándar contra los datos:
>#normal assumption on logit scale, back transform to probability scales
>pred_norm<-predict(glm1,se.fit=T)
>lower<-pred_norm$fit-1.96*pred_norm$se.fit
>upper<-pred_norm$fit+1.96*pred_norm$se.fit
>y<-cbind(lower,pred_norm$fit,upper)
>ypred_norm<-1./(1+exp(-y))*nest_data$n
>with(nest_data,plot(elev,y))
>with(nest_data,matlines(elev,ypred_norm,col="red",lty=c(1,2,2)))
Por último, cabe destacar que glm () proporciona una estadística adicional llamado AIC o el criterio de información de Akaike. AIC es una métrica, basado en la verosimilitud, calculada para un modelo específico dado los datos,
que vamos a utilizar para la comparación entre modelos alternativos más adelante en este laboratorio.
Modelos de regresión múltiple
Modelos de regresión múltiple funcionan exactamente como la regresión lineal y modelos lineales generalizados con predictores individuales, sólo que ahora tenemos más de una variable de predicción. Así que ahora el modelo es de la forma general
Y = b0 + b1 X1 + b2X2 +b3X3...... + bkXk
donde Y es de nuevo la respuesta prevista y X1.... Xk son predictores en cada muestra. Finalmente, b0, b1 ..... bk son los coeficientes del modelo o parámetros que se estiman ya sea utilizando mínimos cuadrados o máxima verosimilitud. Hay muchas formas de que el lineal múltiple modelo puede tomar, al igual que siempre son lineales en los parámetros. Al igual que con los modelos lineales simples, el
respuesta, predictores, o ambos pueden ser transformados antes de ajuste del modelo utilizando logarítmica, logit, sinusoidal, raíz cuadrada, u otras transformaciones. Algunos de los predictores en el modelo puede implicar ya sea en términos de polinomios o de interacción. Por ejemplo
Y = b0 + b1 X+b2 X2
sólo tiene un factor de predicción medido, pero tiene tanto un un término cuadrático lineal y (así, esto sería un ejemplo de regresión polynomial.
Ejemplo:
Aquí tenemos datos de 70 nidos localizados en una sola temporada. El archivo de datos contiene la respuesta (en este caso, 1 = el nido sucedió 0 = error, la elevación en metros, de fecha Julian (días después del 1 de enero) la se encontró el nido, la abundancia oruga (número de orugas / 50 hojas de los árboles), y abundancia de arbustos (número de tallos por 3m X 3m trama centrada en el nido). Podemos empezar por la lectura de los datos y la construcción de un modelo que predice la respuesta en términos de elevación y la fecha juliana. Los datos son aquí y el código de R está aquí
> #MULTIPLE REGRESSION
> nest_multiple<-read.csv("2011_nest_data.csv")
> glm1<-glm(success~elev+jdate,family=binomial(link=logit), data=nest_multiple)
> summary(glm1)
Tenga en cuenta que en este ejemplo no estamos guardando un registro de número de nidos, porque en cada lugar sólo hay 1 nido que tiene éxito (1) o no (0). Por lo tanto, no necesitamos un factor de ponderación que en el ejemplo anterior binomial.
Los resultados resumidos para el GLM son
Call:
glm(formula = success ~ elev + jdate, family = binomial(link = logit))
Deviance Residuals:
Min 1Q Median 3Q Max
-0.9235 -0.5698 -0.3120 -0.2000 2.4537
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.126836 5.270897 0.593 0.5530
elev 0.004478 0.003330 1.345 0.1786
jdate -0.062751 0.028466 -2.204 0.0275 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 49.754 on 69 degrees of freedom
Residual deviance: 42.065 on 67 degrees of freedom
AIC: 48.065
Number of Fisher Scoring iterations: 6
Así, en la escala logit este modelo está prediciendo 1) una intercepción (ambos predictores = 0) de 3,12, un aumento de 0.004 con la elevación, y una disminución de 0.06 con cada día Juliano. Graficar de la respuesta predicha con los datos es un poco complicado para este ejemplo, por 2 razones. En primer lugar, la respuesta es un 1 (éxito) o 0 (fallo) pero realmente predicen la probabilidad de éxito
o el fallo, por lo tanto, no podemos realmente esperar las predicciones del modelo se alineen con los datos. La otra cuestión es que si deseamos mostrar las predicciones en 2 dimensiones (el diagrama de dispersión usual enfoque) tenemos que elegir uno predictor como el eje x y medio sobre el segundo. Así por ejemplo, para trazar contra un rango de la elevación con el promedio de fecha juliana podemos usar un código como éste
> elev0<-seq(900,1400,50)
> jdate<-mean(nest_multiple$jdate)
> new <- data.frame(elev=elev0,jdate=jdate)
> ypred2<-predict(glm1, new=new, se.fit=T)
> #approximate 95% ci on logit scale
> est<-ypred2$fit
> lower<-ypred2$fit-1.96*ypred2$se.fit
> upper<-ypred2$fit+1.96*ypred2$se.fit
> #inverse logit transform
> expit<-function(x){1/(1+exp(-x))}
> ypred<-data.frame(est=expit(est),lower=expit(lower),upper=expit(upper),elev=new$elev)
> with(ypred,plot(elev,est,type="l",ylim=c(0,0.6)))
> with(ypred,matlines(elev,lower,lty=2))
> with(ypred,matlines(elev,upper,lty=2))
Esto muestra un gráfico de la respuesta pronosticada (una probabilidad) y en base normal-(en el logit escala) los límites de confianza, superpuestos sobre la respuesta de datos. Esencialmente podríamos revertir este procedimiento para examinar la variación en el éxito en relación a Julian fecha, un promedio de más altitud.
Evaluar y comparar múltiples modelos
El ejemplo anterior nos conduce rápidamente a un problema (o una oportunidad, si quieres ver de esa manera): tenemos más de un posible modelo que se puede utilizar para describir los datos. Así, por encima de hemos construido un modelo que predice el éxito de nidificación como una función lineal de la elevación y la fecha juliana:
> glm1<-glm(success~elev+jdate,family=binomial(link=logit),data=nest_multiple)
Sin embargo, con la misma facilidad podría haber postulado modelos, puesto que sólo la elevación, fecha juliana
solo, o fecha juliana interactuar con la elevación, predecían el éxito, así que:
> glm2<-glm(success~elev,family=binomial(link=logit),data=nest_multiple)
> glm3<-glm(success~jdate,family=binomial(link=logit),data=nest_multiple)
> glm4<-glm(success ~ elev|jdate, family = binomial(link = logit),data=nest_multiple)
("Elev | jdate" es la abreviatura de "elev + jdate + elev * jdate"
También podríamos haber entretenido a un modelo sin predictores (sólo una intercepción)
> glm5<- glm(formula = success ~ 1, family = binomial(link = logit),data=nest_multiple)
("1" aquí es "uno" no minúscula "el")
Así que ahora tenemos 5 modelos candidatos, y que ni siquiera hemos hablado acerca de la abundancia o de la oruga la densidad de arbustos como posibles predictores en sustitución de o además de la elevación y la fecha juliana. ¿Cómo decidimos qué modelo utilizar? Hay varios criterios que deben tenerse en cuenta en
la construcción de modelos y comparación, así que vamos a ir a través de ellos.
Plausibilidad de los modelos
Soy un firmo defensor de la construcción de modelos que tienen sentido biológico y que tienen tras de sí una hipótesis biológica plausible. Es muy fácil de construir modelos complejos que implica muchos factores, pero muchos de estos modelos es probable que sean sin sentido (o al menos no hay plausibles mecanismo detrás de la modelo). Para volver al ejemplo de éxito de los nidos, tenemos hipótesis plausibles de por qué el éxito de nidificación se prevé que aumente con la elevación (que tiene que ver con calendario de la disponibilidad de alimentos en relación con la fenología del nido) y una disminución con fecha juliana de inicio(Disminución general de la disponibilidad de alimentos y el aumento de los depredadores como avanza la temporada). Podemos también hacer un buen caso para una interacción entre estos factores. Por lo tanto, para cada uno de los modelos
> glm1<-glm(success~elev+jdate,family=binomial(link=logit),data=nest_multiple)
> glm2<-glm(success~elev,family=binomial(link=logit),data=nest_multiple)
> glm3<-glm(success~jdate,family=binomial(link=logit),data=nest_multiple)
> glm4<-glm(success ~ elev|jdate, family = binomial(link = logit),data=nest_multiple)
podemos afirmar un mecanismo plausible. De hecho, el desarrollo de estos modelos debería, y de hecho hizo, se producen en forma inversa: que indica la hipótesis biológica, y luego la conversión a un modelo predictivo. Por otro lado, no podemos imaginar una hipótesis plausible que implica la interacción de la elevación y la fecha juliana en la ausencia del efecto principal de cualquiera de los factores. Así
el modelo
> glm6<- glm(formula = success ~ I(sqrt(elev*jdate)), family = binomial(link = logit),data=nest_multiple)
>
I(sqrt(elev*jdate)) indica que se usa el raiz cuadrado de elev*jdata como la variable
hace que (para nosotros) no tiene sentido. Nosotros no construir este modelo, ya que por una cosa que no tenemos idea cómo interpretar sus estimaciones o predicciones.
Ajuste del model
Suponiendo que nos hemos limitado a un conjunto de modelos plausibles, tenemos una serie de criterios que se pueden utilizar para evaluar los distintos modelos y comparar entre modelos. Ya tenemos visto varias medidas de ajuste del modelo. La mayoría de las medidas estándar de ajuste del modelo de regresión se basan en la suma residual de cuadrados (RSS) también conocida como la desviación, con la idea de que menor RSS se traduce en un mejor ajuste del modelo a los datos. Otras medidas, como el R2, la proporción de la variación respuesta explicada por el modelo, y los errores estándar de estimados de los parámetros todos están directamente relacionados con RSS, con la desviación de máxima verosimilitud -2 log-verosimilitud
en su turno son derivados directemente de la función de verosimilitud. Por lo tanto, una aparente respuesta a la pregunta "¿Cuál es el mejor modelo?" Sería "El que minimiza MSS o maximiza R2 ". Sin embargo, podemos hacer que el flujo lo más pequeño posible (todos los manera de cero), y al mismo tiempo maximizar
R2 (a 1 si se quiere) simplemente añadiendo más términos al modelo lineal. Podemos ilustrar esto con un ejemplo simple de n = 10 observaciones de Y y
X, donde estamos de montaje serie de modelos lineales cada vez más complejos con los términos del polinomio. Los datos están aquí y los comandos en R están aquí
>
Podemas ajustar un serie de modelos más complejos con términos polinomiales.
> #Polynomial regression example of overfitting
> poly_data<-read.csv("polynomial_example.csv")
>
> reg1<-lm(y~x,data=poly_data)
> reg2<-lm(y~poly(x,2),data=poly_data)
> reg3<-lm(y~poly(x,3),data=poly_data)
> reg4<-lm(y~poly(x,4),data=poly_data)
> reg5<-lm(y~poly(x,5),data=poly_data)
> reg6<-lm(y~poly(x,6),data=poly_data)
> reg7<-lm(y~poly(x,7),data=poly_data)
> reg8<-lm(y~poly(x,8),data=poly_data)
> reg9<-lm(y~poly(x,9),data=poly_data)
donde
lm (y ~ x, data = poly_data)
crea el modelo
Y =b0+b1X
reg2<-lm(y~poly(x,2,data=poly_data)) creates
crea
Y =b0+b1X+b2X2
etc.
Supongamos que construimos todos los modelos a un polinomio de orden 9 ª, montarlas a los datos, y utilizar los modelos resultantes para predecir Y en X = 70 (a la derecha en el medio de la gama de valores de X observados).
Podemos resumir los resultados de los modelos de 9 de esta tabla:
model R_squared RSS sigma Loglik n_par Pred se
1 reg1 0.7072946 482.832561 7.768788 -33.57481 2 49.18729 2.811745
2 reg2 0.8005478 329.006647 6.855724 -31.65685 3 52.07049 2.949019
3 reg3 0.8166699 302.412446 7.099442 -31.23542 4 56.86712 7.275362
4 reg4 0.8344578 273.070293 7.390133 -30.72511 5 48.69590 13.477015
5 reg5 0.8832884 192.521767 6.937611 -28.97751 6 47.85469 12.668473
6 reg6 0.8923382 177.593669 7.694016 -28.57395 7 78.22444 62.087778
7 reg7 0.8999276 165.074516 9.085002 -28.20844 8 202.82602 328.226778
8 reg8 0.9956562 7.165257 2.676800 -12.52268 9 409.40665 106.249619
9 reg9 1.0000000 0.000000 NaN Inf 10 755.62975 NaN
Ustedes deben notar varios patrones con estos resultados
A medida que aumentaba la complejidad del modelo, aumenta - todo el camino a 1 para los más complejos modelos.
RSS disminuye y aumenta el logaritmo de verosimilitud con la complejidad, la RSS de ir a cero (Indicando ajuste perfecto) y log-verosimilitud al infinito. Hmm
sigma se hace más pequeño, pero luego se convierten en indefinidos porque para el modelo reg9 se trata dividir RSS por n = 0 -10 grados de libertad.
El error de predicción generalmente aumenta con la complejidad y, finalmente, se convierte en indefinido (el problema anterior).
Lo que está ocurriendo es la disyuntiva clásica en la modelización estadística entre sesgo y la precisión.
Para un conjunto dado de datos tenemos una cantidad fija de "contenido de la información"
A medida que aumentamos la complejidad del modelo mejoramos sesgo de modelo y de mejorar en forma (así, inferior RSS y más alto
Sin embargo, a medida que añadimos parámetros y complejidad (pero no agregue más información a partir de datos) se obtiene estimaciones menos exactas de los parámetros en el modelo
En consecuencia, la fiabilidad de los estimados de los parámetros, y en especial de las predicciones de nuevas observaciones, tenderán a deteriorarse si los modelos se vuelven demasiado complejos
Hay una serie de aproximaciones que se han implementado para optimizar la complejidad en modelos lineales. Vamos a considerar un enfoque más "tradicional", y luego haremos hincapié un enfoque que nos parece superior, basado en la teoría de la información.
Aproximación de las pruebas de hipótesis
El enfoque que se ha defendido es el ajuste de los modelos que están secuencialmente más complejo, con modelos más sencillos anidados en otros más complejos, y luego utilizando la prueba de hipótesis para comparar modelos más simples (null) vs más complejas (alternativos), por ejemplo a través de pruebas de coeficiente de riesgo.
Por ejemplo, el modelo de
Y =b0+b1X
se puede considerar como una simplificación del modelo (H 1)
Y =b0+b1X+b2X2
bajo la hipótesis nula.
H0:b2 = 0
Una prueba de razón verosimilitud se llevaría a cabo usando las verosimilitudes para cada modelo.
Esta estadística sigue una distribución chi-cuadrado bajo H 0, con grados de libertad igual a la
diferencia en el número de parámetros entre H 1 y H 0. Para el ejemplo anterior, lo haríamos
calcular esto como
> #LRT
> T<--2*(logLik(reg1)[1]-logLik(reg2)[1])
> df<-summary(reg2)$df[1]-summary(reg1)$df[1]
> p<-pchisq(T,df,lower.tail=F)
> T
[1] 3.83592
> df
[1] 1
> p
[1] 0.05016546
Esto es justo en el límite para una prueba con el error de tipo I = 0,05, por lo que tal vez le eligió
mantener el modelo más complejo. Entonces se puede proceder de forma secuencial, de dejar de fumar cuando el modelo nulo (ahora reg2) ya no es rechazada (por ejemplo, en comparación con REG1). Así, la repetición de estos cálculos, pero ahora con reg2 como nula y reg3 como la alternativa que tenemos
> T<--2*(logLik(reg2)[1]-logLik(reg3)[1])
> df<-summary(reg3)$df[1]-summary(reg2)$df[1]
> p<-pchisq(T,df,lower.tail=F)
> T
[1] 0.4623851
> df
[1] 1
> p
[1] 0.4965114
así que podríamos suponer que parar en reg2, al no rechazar ese modelo frente a la siguiente más
alternativa complejo.
Aunque las pruebas de hipótesis nula puede trabajar para la selección de modelos, hay una serie de cuestiones, no el menor de ellos es lo que pide es la prueba de hecho (más simple al más complejo, a la inversa, o algo más). Prueba de hipótesis nula puede llevar a overfitting de modelos. También es obviamente limitada a los casos en que el modelo más simple está anidado como hipótesis nula en el más complejo modelo. Sin embargo, hay muchos casos en que el modelo más simple no está anidado en el más uno complejo. Por ejemplo, podemos especificar
H1: Y =b0+b1X+b2X2
H2: Y =b0+b2X3
Ningun de estos se pueden formar como nula del otro, aunque ambas son hipótesis nulas cuando
en comparación con
H3: Y =b0+b1X+b2X2 +b3X3
En su lugar defendemos el enfoque desarrollado por Akaike (1973) y se describen en detalle en el
excelente libro de Burnham y Anderson (2002). El enfoque utiliza la teoría de la información a
desarrollar una estadística conocida como Criterio de Información de Akaike (AIC). AIC se calcula a partir de la log-verosimilitud y el número de parámetros del modelo por la fórmula
AIC = - 2ln log L + 2 k
donde k es el número de parámetros estimados en el modelo (por ejemplo, k = 2 para una regresión simple: b0 y b1. Por lo general, vamos a querer ajustar AIC para el tamaño de la muestra (en especial por corrección para los pequeños muestras). El ajuste será especialmente importante como k se hace grande con respecto a n (como en el ejemplo de regresión polinomio). Volviendo al ejemplo de que podemos calcular AIC c para cada modelo. El valor real de AIC no es importante, más bien es el valor de un modelo de valor de AIC con respecto a la del modelo con la AIC más pequeño. De este modo se forma la diferencia entre los modelos de AIC (AIC c ).
D AIC =AIC - min(AIC)
donde min (AIC) es el valor mínimo AIC sobre el conjunto de los modelos considerados. Finalmente,
D AIC se puede utilizar para calcular un "peso de la información" para el modelo, que representa en esencia la forma tanto que modelo debe contribuir a la inferencia, dados los datos.
A continuación calculamos los pesos ajustados AIC para el ejemplo de polinomio y añadimos
esta información a la tabla. Esto lo hacemos sólo para los primeros 7 modelos, ya ajustado AIC no puede ser calculado para los últimos modelos de 2, debido a la pérdida de grados de libertad.
modelo lik n_par aic_c delta_aic peso
1 reg1 -39.35345 86.42119 0.000000 2 8.229101e-01
AIC nk ls e log (2) 2 = s +
Hemos aplicado esta aproximación para comparar los modelos del ejemplo de regresión polinómica, con los resultados siguintes
> models2
model lik n_par aic_ls_c delta_aic wt
1 reg1 -39.35345 2 58.27385 0.000000 9.035092e-01
2 reg2 -38.80857 3 62.80512 4.531264 9.375215e-02
3 reg3 -38.57738 4 69.88424 11.610385 2.721257e-03
4 reg4 -38.22156 5 79.99582 21.721966 1.734074e-05
5 reg5 -36.86750 6 94.51913 36.245278 1.217224e-08
6 reg6 -36.79116 7 127.24328 68.969423 9.537018e-16
7 reg7 -36.78487 8 221.28534 163.011490 3.617739e-36
Los resultados muestran ordenamiento similar de mínimos cuadrados basado AIC con probabilidad AIC, pero con considerablemente más peso en el modelo reg1. Recomendamos el uso de AIC probabilidad basada en cuando posible, pero si no hay un modelo de probabilidad está disponible, por mínimos cuadrados AIC es un razonablemente alternativa.
Vea el código de secuencia de comandos R para los detalles de cómo se realizan los cálculos anteriores.
AIC es una medida del contenido de información de los datos frente a la complejidad del modelo, y tiene ha demostrado resultar en la selección del modelo óptimo, sobre todo cuando el objetivo es la predicción fuera de los datos del modelo. Más tarde, en el laboratorio se considera una alternativa (o en realidad un suplemento) para AIC, con base en el ordenador intensiva de re-muestreo de los datos.
Inferencia y el modelo multi promedio de modelo
En lo anterior se hizo hincapié en la comparación entre los modelos, con la idea de la selección de solo un modelo (o tal vez 1 o 2 modelos) en la que a la inferencia base. Sin embargo, a menudo es el caso que un parámetro o predicción se produce en varios modelo. Si nuestro objetivo principal es la estimación o predicción, entonces vamos a preocuparse más por la estimación adecuada de ese parámetro (o predicción) que sobre lo que el modelo o modelos que se basa en. Supongamos que θ es un parámetro o una predicción que estamos interesados en la estimación, y que se dispone de estimaciones con arreglo a varios candidatos modelos. Burnham y Anderson (2002) describen cómo podemos combinar los estimados sobre el modelos alternativos utilizando los pesos de Akaike.
La inferencia también requiere una estimado de la precisión del estimador. Debajo de cada modelo de candidato hay un error estándar estimado (según sea la desviación o el máximo de mínimos cuadrados
probabilidad),. Sin embargo, este estimdor, por definición, se supone que el modelo de la estimación es el verdadero modelo. Para dar cuenta de la incertidumbre en la estimación de parámetros debido a la incertidumbre del modelo, Burnham y Anderson (2002) cálcularon un "error no condicionado estándar". En general, el uso del error estándar incondicional producirá mayores estimados de varianza y intervalos de confianza más amplios en consecuencia, lo que es apropiado dado que la inferencia no descansa en la creencia en un solo modelo.
Ejemplo
Volvemos al ejemplo de éxito de los nidos y consideramos 3 modelos alternativos: el éxito de nidificación como una función de la elevación y la fecha juliana, y los modelos con cada uno de estos factores solamente. También consideradamos un modelo con la elevación, fecha juliana, y la interacción de estos factores, pero esto modelo no produjeron estimaciones razonables y por lo tanto fue excluido. Estamos interesados específicamente en la predicción de la respuesta (éxito de los nidos en la escala logit y luego de vuelta transformado de probabilidad) en la cota 1200m y Julian Fecha = 176 (cerca de los medios para estos valores en el estudio). Corrimos todos los 3 de estos modelos en glm () y calculamos predicciones, errores típicos, y estadísticas de la AIC. Los valores se utilizan entonces para producir una predicción de modelo promediado, incondicional SE, y 95% intervalos de confianza normal en la escala logit, que eran la parte de atrás-transformado para producir estimados sobre la escala de probabilidad.
> #predictions
> p1<-predict(glm1,list(elev=1200,jdate=176),se.fit=TRUE)
> p2<-predict(glm2,list(elev=1200,jdate=176),se.fit=TRUE)
> p3<-predict(glm3,list(elev=1200,jdate=176),se.fit=TRUE)
> #p4<-predict(glm4,list(elev=1200,jdate=176),se.fit=TRUE)
> pred<-c(p1$fit,p2$fit,p3$fit)
> se_pred<-c(p1$se.fit,p2$se.fit,p3$se.fit)
> #model labels
> mods<-c("elev+jdate","elev","jdate")
> npar<-c(summary(glm1)$df[1],summary(glm2)$df[1],summary(glm3)$df[1])
> n<-nrow(nest_multiple)
> aic<-c(AIC(glm1),AIC(glm2),AIC(glm3))
> aic_c=aic+2*npar*(npar+1)/(n-npar-1)
> delta<-aic_c-min(aic_c)
> Z<-exp(-0.5*delta)
> wt<-Z/sum(Z)
>model_avg_pred<-sum(wt*pred)
>uncond_se<-sum(wt*sqrt(se_pred^2+(pred-model_avg_pred)^2))
> model_avg_pred
[1] -2.447551
> uncond_se
[1] 0.5378216
> print("predicted probabilty")
[1] "predicted probabilty"
> 1./(1+exp(-model_avg_pred))
[1] 0.07961781
> lower<-model_avg_pred-1.96*uncond_se
> upper<-model_avg_pred+1.96*uncond_se
> print("CI")
[1] "CI"
> c(1./(1+exp(-lower)), 1./(1+exp(-upper)))
[1] 0.02926443 0.19886219
Los datos para este ejemplo
Los comandos de R para el ejemplo
Siguiente: Escribir funciones definido por el usuario