Lab 7. Modelación no lineal
En este laboratorio, volvemos a analizar los datos de Wage considerados en los ejemplos a lo largo de este capítulo, con el fin de ilustrar el hecho de que muchos de los complejos procedimientos de ajuste no lineales discutidos pueden ser fácilmente implementados en R. Comenzamos cargando la librería ISLR, que contiene los datos.
Link de Rpubs: https://rpubs.com/arturolozanoo/1194138
library (ISLR)
attach (Wage)
summary (Wage)## year age maritl race
## Min. :2003 Min. :18.00 1. Never Married: 648 1. White:2480
## 1st Qu.:2004 1st Qu.:33.75 2. Married :2074 2. Black: 293
## Median :2006 Median :42.00 3. Widowed : 19 3. Asian: 190
## Mean :2006 Mean :42.41 4. Divorced : 204 4. Other: 37
## 3rd Qu.:2008 3rd Qu.:51.00 5. Separated : 55
## Max. :2009 Max. :80.00
##
## education region jobclass
## 1. < HS Grad :268 2. Middle Atlantic :3000 1. Industrial :1544
## 2. HS Grad :971 1. New England : 0 2. Information:1456
## 3. Some College :650 3. East North Central: 0
## 4. College Grad :685 4. West North Central: 0
## 5. Advanced Degree:426 5. South Atlantic : 0
## 6. East South Central: 0
## (Other) : 0
## health health_ins logwage wage
## 1. <=Good : 858 1. Yes:2083 Min. :3.000 Min. : 20.09
## 2. >=Very Good:2142 2. No : 917 1st Qu.:4.447 1st Qu.: 85.38
## Median :4.653 Median :104.92
## Mean :4.654 Mean :111.70
## 3rd Qu.:4.857 3rd Qu.:128.68
## Max. :5.763 Max. :318.34
## 1. Regresión Polinómica y Funciones de paso
Para este ejemplo ajustaremos a un modelo lineal, usando la función lm(), para predecir el salario usando un polinomio de cuarto grado en age: poly(age,4). El comando poly() nos permite evitar tener que escribir una fórmula larga. El comando es el siguiente:
fit<-lm(wage~poly(age,4), data=Wage)
summary(fit)##
## Call:
## lm(formula = wage ~ poly(age, 4), data = Wage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -98.707 -24.626 -4.993 15.217 203.693
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 111.7036 0.7287 153.283 < 2e-16 ***
## poly(age, 4)1 447.0679 39.9148 11.201 < 2e-16 ***
## poly(age, 4)2 -478.3158 39.9148 -11.983 < 2e-16 ***
## poly(age, 4)3 125.5217 39.9148 3.145 0.00168 **
## poly(age, 4)4 -77.9112 39.9148 -1.952 0.05104 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 39.91 on 2995 degrees of freedom
## Multiple R-squared: 0.08626, Adjusted R-squared: 0.08504
## F-statistic: 70.69 on 4 and 2995 DF, p-value: < 2.2e-16La función devuelve una matriz cuyas columnas son una base de polinomios ortogonales, lo que esencialmente significa que cada columna es una combinación lineal de las variables \(age\), \(age^{2}\), \(age^{3}\) y \(age^{4}\).
Ahora creamos una cuadrícula de valores para age en la que queremos predicciones, y luego llamamos a la función genérica predict(), especificando que también queremos errores estándar.
agelims<-range(age)
agelims## [1] 18 80age.grid<-seq (from=agelims [1], to=agelims [2])
age.grid## [1] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
## [26] 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
## [51] 68 69 70 71 72 73 74 75 76 77 78 79 80preds<-predict(fit, newdata=list(age=age.grid),se=TRUE)
preds## $fit
## 1 2 3 4 5 6 7 8
## 51.93145 58.49674 64.57188 70.18273 75.35440 80.11122 84.47676 88.47380
## 9 10 11 12 13 14 15 16
## 92.12437 95.44973 98.47038 101.20601 103.67560 105.89731 107.88856 109.66599
## 17 18 19 20 21 22 23 24
## 111.24548 112.64214 113.87029 114.94351 115.87459 116.67557 117.35770 117.93148
## 25 26 27 28 29 30 31 32
## 118.40662 118.79210 119.09608 119.32598 119.48846 119.58939 119.63388 119.62627
## 33 34 35 36 37 38 39 40
## 119.57013 119.46827 119.32271 119.13473 118.90481 118.63269 118.31733 117.95691
## 41 42 43 44 45 46 47 48
## 117.54885 117.08980 116.57566 116.00152 115.36174 114.64988 113.85877 112.98043
## 49 50 51 52 53 54 55 56
## 112.00613 110.92638 109.73090 108.40865 106.94784 105.33588 103.55943 101.60437
## 57 58 59 60 61 62 63
## 99.45583 97.09815 94.51491 91.68893 88.60223 85.23611 81.57105
##
## $se.fit
## 1 2 3 4 5 6 7 8
## 5.298268 4.370763 3.592101 2.955813 2.455588 2.083260 1.825676 1.662329
## 9 10 11 12 13 14 15 16
## 1.566864 1.512533 1.477572 1.447355 1.413769 1.373549 1.326690 1.275230
## 17 18 19 20 21 22 23 24
## 1.222382 1.171865 1.127324 1.091786 1.067171 1.053981 1.051263 1.056899
## 25 26 27 28 29 30 31 32
## 1.068100 1.081939 1.095809 1.107729 1.116537 1.121992 1.124827 1.126747
## 33 34 35 36 37 38 39 40
## 1.130363 1.139015 1.156439 1.186260 1.231415 1.293640 1.373218 1.469069
## 41 42 43 44 45 46 47 48
## 1.579082 1.700568 1.830715 1.966988 2.107496 2.251362 2.399124 2.553187
## 49 50 51 52 53 54 55 56
## 2.718307 2.902040 3.115018 3.370894 3.685814 4.077417 4.563601 5.161438
## 57 58 59 60 61 62 63
## 5.886530 6.752911 7.773327 8.959688 10.323512 11.876272 13.629642
##
## $df
## [1] 2995
##
## $residual.scale
## [1] 39.91479se.bands<-cbind(preds$fit + 2*preds$se.fit, preds$fit - 2*preds$se.fit)Finalmente, ploteamos los datos y añadimos el ajuste del polinomio de grado 4.
par(mfrow =c(1,1), mar=c(4.5 ,4.5 ,1 ,1), oma=c(0,0,4,0))
plot(age, wage, xlim=agelims, cex =.5, col="darkgrey")
title ("Polinomio de 4to grado", outer =T)
lines(age.grid, preds$fit, lwd=2, col=" blue")
matlines (age.grid, se.bands, lwd=1, col=" blue", lty =3)
Aquí los argumentos mar y oma de par() nos permiten controlar los márgenes del plot, y la función title() crea un título de figura que abarca ambos ploteos.
Al realizar una regresión polinómica debemos decidir el grado del polinomio a utilizar. Una forma de hacerlo es usando pruebas de hipótesis. En este caso ajustaremos modelos que van desde el polinomio lineal hasta un grado 5 y tratamos de determinar el modelo más simple que sea suficiente para explicar la relación entre wage y age. Utilizamos la función anova(), que realiza un análisis de varianza (ANOVA por sus siglas en inglés, utilizando una prueba F) para probar la hipótesis nula de que un modelo M1 es suficiente para explicar los datos frente a la hipótesis alternativa de que se requiere un modelo M2 más complejo. Para poder utilizar la función anova(), M1 y M2 deben ser modelos anidados: los predictores en M1 deben ser un subconjunto de los predictores en M2. En este caso, ajustamos cinco modelos diferentes y comparamos secuencialmente el modelo más simple con el más complejo.
fit.1 <- lm(wage~age,data=Wage)
fit.2 <- lm(wage~poly(age,2),data=Wage)
fit.3 <- lm(wage~poly(age,3),data=Wage)
fit.4 <- lm(wage~poly(age,4),data=Wage)
fit.5 <- lm(wage~poly(age,5),data=Wage)
anova(fit.1, fit.2, fit.3, fit.4, fit.5)## Analysis of Variance Table
##
## Model 1: wage ~ age
## Model 2: wage ~ poly(age, 2)
## Model 3: wage ~ poly(age, 3)
## Model 4: wage ~ poly(age, 4)
## Model 5: wage ~ poly(age, 5)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2998 5022216
## 2 2997 4793430 1 228786 143.5931 < 2.2e-16 ***
## 3 2996 4777674 1 15756 9.8888 0.001679 **
## 4 2995 4771604 1 6070 3.8098 0.051046 .
## 5 2994 4770322 1 1283 0.8050 0.369682
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El p-value que compara el Modelo 1 lineal con el Modelo 2 cuadrático es esencialmente cero (<\(10^{-16}\)), lo que indica que un ajuste lineal no es suficiente. De manera similar, el p-value que compara el Modelo 2 cuadrático con el Modelo 3 cúbico es muy bajo (0,001679), por lo que el ajuste cuadrático también es insuficiente. El p-value que compara los polinomios cúbicos y de grado 4, Modelo 3 y Modelo 4, es aproximadamente el 5% mientras que el polinomio de grado 5 Modelo 5 parece innecesario porque su p-value es 0.37. Por lo tanto, un polinomio cúbico o de grado 4 parece proporcionar un ajuste razonable a los datos, pero no se justifican los modelos de orden inferior o superior.
A continuación consideramos la tarea de predecir si un individuo gana más de 250.000 dólares al año. Procedemos de la misma manera que antes, excepto que primero creamos el vector de respuesta apropiado, y luego aplicamos la función glm() usando family=“binomial” para ajustar un modelo de regresión logística polinómica.
fit<-glm(I(wage >100)~poly(age,4), data=Wage, family=binomial)Obsérvese que volvemos a utilizar el I() a forma de envoltorio, para crear esta variable de respuesta binaria sobre la marcha. La expresión wage>250 se evalúa a una variable lógica que contiene TRUEs y FALSEs, que glm() coacciona a binario poniendo los TRUEs a 1 y los FALSEs a 0.
Una vez más, hacemos predicciones usando la función predict().
preds<-predict(fit, newdata=list(age=age.grid), se=T)Sin embargo, el cálculo de los intervalos de confianza está ligeramente más complicado que en el caso de la regresión lineal. El tipo de predicción por defecto para un modelo glm() es type=“link”, que es lo que usamos aquí. Esto significa que obtenemos predicciones para el logit: es decir, hemos ajustado un modelo de la forma
\[log\frac{Pr(Y=1|X)}{1-Pr(Y=1|X)}=X\beta,\]
y las predicciones que se han hecho son de la forma \(X\hat{\beta}\) que se ha hecho. Los errores estándar dados también son de esta forma. Para obtener los intervalos de confianza para Pr(\(Y\) = 1|\(X\)), usamos la transformación
\[Pr(Y=1|X)=\frac{exp(X\beta)}{1+exp(X\beta)}.\]
pfit<-exp(preds$fit)/(1 + exp( preds$fit ))
se.bands.logit <- cbind(preds$fit + 2*preds$se.fit, preds$fit - 2*preds$se.fit)
se.bands <- exp(se.bands.logit)/(1+ exp(se.bands.logit))Finalmente, se grafica de la siguiente manera:
plot(age, I(wage>100), xlim=agelims, type="n", ylim=c(0,1))
points(jitter(age), I((wage>250)/5), cex=0.5, pch="|", col="darkgrey")
lines(age.grid, pfit, lwd=2, col="blue")
matlines(age.grid, se.bands, lwd=1, col="blue", lty =3)
Hemos dibujado los valores de edad correspondientes a las observaciones con valores de wage superiores a 250 como marcas grises en la parte superior del gráfico, y los que tienen valores de wage inferiores a 250 se muestran como marcas grises en la parte inferior del gráfico. Hemos utilizado la función jitter() para hacer oscilar un poco los valores de wage de modo que las observaciones con el mismo valor de wage no se cubran entre sí. Esto se denomina a menudo un rug plot.
Para ajustar una función de paso, usamos la función cut().
table(cut(age,4))##
## (17.9,33.5] (33.5,49] (49,64.5] (64.5,80.1]
## 750 1399 779 72fit<-lm(wage~cut(age,4), data=Wage)
coef(summary(fit))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94.158392 1.476069 63.789970 0.000000e+00
## cut(age, 4)(33.5,49] 24.053491 1.829431 13.148074 1.982315e-38
## cut(age, 4)(49,64.5] 23.664559 2.067958 11.443444 1.040750e-29
## cut(age, 4)(64.5,80.1] 7.640592 4.987424 1.531972 1.256350e-01Aquí cut() automáticamente escogió los puntos de corte a los 33.5, 49 y 64.5 años de edad. También podríamos haber especificado nuestros propios puntos de corte directamente usando la opción de breaks. La función cut() devuelve una variable categórica ordenada; la función lm() crea entonces un conjunto de variables ficticias para su uso en la regresión. La categoría age<33.5 se deja fuera, por lo que el coeficiente del intercepto de $94,160 dólares puede interpretarse como el salario medio para los menores de 33.5 años, y los otros coeficientes pueden interpretarse como el salario adicional medio para los de los otros grupos de edad. Podemos producir predicciones y gráficos como lo hicimos en el caso del ajuste polinómico.
2. Splines
Para ajustar las regresiones splines en R, usamos la libreria de splines. las regresiones splines pueden ser ajustadas construyendo una matriz apropiada de funciones base. La función bs() genera toda la matriz de funciones base para splines con el conjunto de nodos especificado. Por defecto, se producen splines cúbicos. El ajuste de wage a age mediante una regresión spline es sencillo:
library(splines)
fit<-lm(wage~bs(age, knots=c(25, 40, 60)), data=Wage)
pred<-predict(fit, newdata=list(age=age.grid), se=T)
plot(age, wage, col="gray")
lines(age.grid, pred$fit ,lwd =2)
lines(age.grid, pred$fit + 2*pred$se, lty="dashed")
lines(age.grid, pred$fit - 2*pred$se, lty="dashed")
summary(fit)##
## Call:
## lm(formula = wage ~ bs(age, knots = c(25, 40, 60)), data = Wage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -98.832 -24.537 -5.049 15.209 203.207
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 60.494 9.460 6.394 1.86e-10 ***
## bs(age, knots = c(25, 40, 60))1 3.980 12.538 0.317 0.750899
## bs(age, knots = c(25, 40, 60))2 44.631 9.626 4.636 3.70e-06 ***
## bs(age, knots = c(25, 40, 60))3 62.839 10.755 5.843 5.69e-09 ***
## bs(age, knots = c(25, 40, 60))4 55.991 10.706 5.230 1.81e-07 ***
## bs(age, knots = c(25, 40, 60))5 50.688 14.402 3.520 0.000439 ***
## bs(age, knots = c(25, 40, 60))6 16.606 19.126 0.868 0.385338
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 39.92 on 2993 degrees of freedom
## Multiple R-squared: 0.08642, Adjusted R-squared: 0.08459
## F-statistic: 47.19 on 6 and 2993 DF, p-value: < 2.2e-16Aquí tenemos nodos preespecificados a las edades de 25, 40 y 60 años. Esto produce un spline con seis funciones básicas. (Recordemos que un spline cúbico con tres nodos tiene siete grados de libertad; estos grados de libertad son usados por un intercepto, más seis funciones básicas). También podríamos usar la opción df para producir una spline con nodos en cuantiles uniformes de los datos.
dim(bs(age, knots=c(25, 40, 60)))## [1] 3000 6dim(bs(age, df=6))## [1] 3000 6attr(bs(age, df=6), "knots")## [1] 33.75 42.00 51.00En este caso R elige nodos a las edades de 33.8, 42.0 y 51.0, que corresponden a los percentiles 25, 50 y 75 de la age. La función bs() también tiene un argumento de degree, por lo que podemos ajustar splines de cualquier grado, en lugar del grado por defecto de 3 (que produce un spline cúbico).
Para poder ajustar un spline natural, usamos la función ns(). Aquí ajustamos un spline natural con cuatro grados de libertad.
fit2<-lm(wage~ns(age, df=2), data=Wage)
pred2<-predict(fit2, newdata=list(age=age.grid), se=T)
plot(age, wage, col="gray")
lines(age.grid, pred2$fit, col="red", lwd =2)
Al igual que con la función bs(), podríamos en cambio especificar los knots directamente usando la opción de nodos.
Para ajustar un spline de alisado, usamos la función smooth.spline(). La figura 7.8 fue producida con el siguiente código:
plot(age, wage, xlim=agelims, cex =.5, col ="darkgrey")
title("Smoothing Spline")
fit<-smooth.spline(age, wage, df =16)
fit2<-smooth.spline(age, wage, cv=TRUE)## Warning in smooth.spline(age, wage, cv = TRUE): cross-validation with
## non-unique 'x' values seems doubtfulfit2$df## [1] 6.794596lines(fit, col ="red", lwd =2)
lines(fit2, col ="blue", lwd =2)
legend("topright", legend=c("16 DF", "6.8 DF"), col=c("red"," blue"), lty =1, lwd =2, cex =.8)
Fíjense que en la primera llamada a smooth.spline(), especificamos df=16. La función entonces determina qué valor de \(\lambda\) lleva a 16 grados de libertad. En la segunda llamada a smooth.spline(), seleccionamos el nivel de suavidad por validación cruzada; esto resulta en un valor de \(\lambda\) que da 6,8 grados de libertad.
3. GAMs (Generalized Additive Models)
Ahora ajustamos un GAM para predecir el wage usando las funciones spline naturales de year y age, tratando education como un predictor cualitativo. Dado que esto es sólo un gran modelo de regresión lineal que utiliza una elección apropiada de funciones base, podemos simplemente hacer esto utilizando la función lm().
gam1<-lm(wage~ns(year, 4) + ns(age, 5) + education, data=Wage)Ahora ajustamos el modelo \(wage=\beta_0+f_1(year)+f_2(age)+f_3(education)+\epsilon\) usando splines suavizados en lugar de splines naturales. Para poder ajustar tipos más generales de GAM, utilizando splines suavizados u otros componentes que no pueden ser expresados en términos de funciones base y luego ajustar usando la regresión de mínimos cuadrados, necesitaremos usar la librería gam en R.
La función s(), que forma parte de la libraría gam, se utiliza para indicar que nos gustaría utilizar un spline suavizado. Especificamos que la función de year debe tener 4 grados de libertad, y que la función de age tendrá 5 grados de libertad. Como education es cualitativa, la dejamos como está, y se convierte en cuatro variables dummy. Utilizamos la función gam() para ajustar un GAM usando estos componentes. Todos los términos en el modelo se ajustan simultáneamente, tomándose en cuenta unos a otros para explicar la respuesta.
library(gam)
gam.m3<-gam(wage~s(year, 4) + s(age, 5) + education, data=Wage)
summary(gam1)##
## Call:
## lm(formula = wage ~ ns(year, 4) + ns(age, 5) + education, data = Wage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -120.513 -19.608 -3.583 14.112 214.535
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 46.949 4.704 9.980 < 2e-16 ***
## ns(year, 4)1 8.625 3.466 2.488 0.01289 *
## ns(year, 4)2 3.762 2.959 1.271 0.20369
## ns(year, 4)3 8.127 4.211 1.930 0.05375 .
## ns(year, 4)4 6.806 2.397 2.840 0.00455 **
## ns(age, 5)1 45.170 4.193 10.771 < 2e-16 ***
## ns(age, 5)2 38.450 5.076 7.575 4.78e-14 ***
## ns(age, 5)3 34.239 4.383 7.813 7.69e-15 ***
## ns(age, 5)4 48.678 10.572 4.605 4.31e-06 ***
## ns(age, 5)5 6.557 8.367 0.784 0.43328
## education2. HS Grad 10.983 2.430 4.520 6.43e-06 ***
## education3. Some College 23.473 2.562 9.163 < 2e-16 ***
## education4. College Grad 38.314 2.547 15.042 < 2e-16 ***
## education5. Advanced Degree 62.554 2.761 22.654 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 35.16 on 2986 degrees of freedom
## Multiple R-squared: 0.293, Adjusted R-squared: 0.2899
## F-statistic: 95.2 on 13 and 2986 DF, p-value: < 2.2e-16Para graficar simplemente llamamos a la función plot():
par(mfrow =c(1,3))
plot(gam.m3, se=TRUE ,col ="blue ")
La función genérica plot() reconoce que gam.m3 es un objeto de la clase gam, e invoca el método plot.Gam() apropiado. Convenientemente, aunque gam1 no es de la clase gam sino de la clase lm, podemos seguir usando plot.Gam() en él.
plot.Gam(gam1, se=TRUE, col ="red")
Noten que tuvimos que usar plot.Gam() en lugar de la función genérica plot().
En estos gráficos, la función de year parece bastante lineal. Podemos realizar una serie de pruebas de ANOVA para determinar cuál de estos tres modelos es el mejor: un GAM que excluye year (M1), un GAM que utiliza una función lineal de year (M2), o un GAM que utiliza una función de spline de year (M3).
gam.m1<-gam(wage~s(age, 5) + education, data=Wage)
gam.m2<-gam(wage~year+s(age, 5) + education, data=Wage)
gam.m3<-gam(wage~s(year, 4) + s(age, 5) + education, data=Wage)
anova(gam.m1, gam.m2, gam.m3, test="F")## Analysis of Deviance Table
##
## Model 1: wage ~ s(age, 5) + education
## Model 2: wage ~ year + s(age, 5) + education
## Model 3: wage ~ s(year, 4) + s(age, 5) + education
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 2990 3711731
## 2 2989 3693842 1 17889.2 14.4771 0.0001447 ***
## 3 2986 3689770 3 4071.1 1.0982 0.3485661
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Encontramos que hay pruebas convincentes de que un GAM con una función lineal de year es mejor que un GAM que no incluye year en absoluto ( p-value=0.00014 ). Sin embargo, no hay evidencia de que se necesite una función no lineal de year ( p-value=0.349 ). En otras palabras, en base a los resultados de este ANOVA, se prefiere M2.
La función summary() produce un resumen del ajuste del GAM.
summary(gam.m2)##
## Call: gam(formula = wage ~ year + s(age, 5) + education, data = Wage)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -119.959 -19.647 -3.199 13.969 213.562
##
## (Dispersion Parameter for gaussian family taken to be 1235.812)
##
## Null Deviance: 5222086 on 2999 degrees of freedom
## Residual Deviance: 3693842 on 2989 degrees of freedom
## AIC: 29885.06
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## year 1 27154 27154 21.973 2.89e-06 ***
## s(age, 5) 1 194535 194535 157.415 < 2.2e-16 ***
## education 4 1069081 267270 216.271 < 2.2e-16 ***
## Residuals 2989 3693842 1236
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## year
## s(age, 5) 4 32.46 < 2.2e-16 ***
## education
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Los p-values de year y age corresponden a una hipótesis nula de una relación lineal frente a la alternativa de una relación no lineal. El elevado p-value para year refuerza nuestra conclusión del test ANOVA de que una función lineal es adecuada para este término. Sin embargo, hay una evidencia muy clara de que se requiere un término no lineal para la age.
Para ajustar una regresión logística GAM, utilizamos una vez más la función I() en la construcción de la variable de respuesta binaria, y establecemos family=binomial.
gam.lr<-gam(I(wage >100)~year+s(age, df=5) + education, family =binomial, data=Wage)
par(mfrow=c(1,3))
plot(gam.lr, se=T, col="green")
Es fácil ver que no hay personas con altos ingresos en la categoría \(<HS\):
table(education, I(wage>250))##
## education FALSE TRUE
## 1. < HS Grad 268 0
## 2. HS Grad 966 5
## 3. Some College 643 7
## 4. College Grad 663 22
## 5. Advanced Degree 381 454. Resuelve
A continuación nos daremos a la tarea de predecir si un individuo gana más de 125,000 dólares al año. Primero creamos el vector de respuesta apropiado, y luego aplicamos la función glm() usando family=“binomial” para ajustar un modelo de regresión logística polinómica.
Realice predicciones usando la función predict().
Calcule los intervalos de confianza
Finalmente, grafica: