Regresiones polinomiales

Cargue de datos

Lo primero que hacemos cada una de las clases es revisar el cargue de nuestros datos;

data <- read.csv("Clase 8.csv", sep=";")
str(data)
## 'data.frame':    352 obs. of  10 variables:
##  $ Marca     : chr  "A Speyside Distillery" "Ainslie's" "Angel's Envy" "Asyla" ...
##  $ Pais      : chr  "Scotland" "Scotland" "United States" "Scotland" ...
##  $ Whiskies  : int  11 32 23 11 12 364 11 12 1071 62 ...
##  $ Votos     : int  10 12 13 7 4 249 11 3 908 26 ...
##  $ Puntuacion: num  86.7 82.3 83.4 81.8 81 ...
##  $ Categoria : chr  "C" "D" "D" "D" ...
##  $ Sabor     : int  8 5 3 3 9 6 3 7 8 7 ...
##  $ Color     : int  6 6 6 2 10 2 7 9 5 9 ...
##  $ Maltas    : int  8 8 1 8 8 1 9 1 1 9 ...
##  $ Compra    : int  0 0 1 0 0 1 0 1 1 0 ...

El día de hoy seguiremos trabajando con la base de datos de las diferentes características de un grupo de whiskies de diferentes países del mundo. La variable que le interesa estimar es el puntaje del whisky. El link para acceder a los datos es el siguiente: https://www.dropbox.com/s/pt6g6luafdyrx9k/Clase%208.csv?dl=0.

Paquete CARET

El paquete caret (classification and regression training, Kuhn (2016)) incluye una serie de funciones que facilitan el uso de decenas de métodos complejos de clasificación y regresión. Utilizar este paquete en lugar de las funciones originales de los métodos presenta dos ventajas:

  • Permite utilizar un código unificado para aplicar reglas de clasificación muy distintas, implementadas en diferentes paquetes.

*Es más fácil poner en práctica algunos procedimientos usuales en problemas de clasificación. Por ejemplo, hay funciones específicas para dividir la muestra en datos de entrenamiento y datos de test o para ajustar parámetros mediante validación cruzada.

Con este paquete nos centraremos en el uso para problemas de regresión y clasificación limitadas, pero son más de 200 modelos diferentes que tendrán a su disposición.

Para ello lo primero que haremos será realizar una partición en los datos en dos grupos, train y test. Esto nos permite entrenar el modelo y, además, probar si funciona. La proporción no es exacta y su elección no suele ser aleatoria, por lo que se aconseja una proporción de 80/20. La función que usaremos es createDataPartition():

library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
trainIndex <- createDataPartition(data$Puntuacion, p=0.8, list=F)

Ahora que tenemos los índices podemos crear nuestros nuevos data frames al usar este índice para hacer la división:

train <- data[trainIndex,]
test <- data[-trainIndex,]

Ahora deberemos definir el resampling para la evaluación del modelo. La finalidad última de un modelo es predecir la variable respuesta en observaciones futuras o en observaciones que el modelo no ha “visto” antes. Para conseguir una estimación más certera, y antes de recurrir al conjunto de test, se pueden emplear estrategias de validación basadas en resampling. caret incorpora los métodos boot, boot632, optimism_boot, boot_all, cv, repeatedcv, LOOCV, LGOCV. Cada uno funciona internamente de forma distinta, pero todos ellos se basan en la idea de ajustar y evaluar el modelo múltiples veces con distintos subconjuntos creados a partir de los datos de entrenamiento, obteniendo en cada repetición una estimación del error. El promedio de todas las estimaciones tiende a converger en el valor real del error de test.

Para especificar el tipo de validación, así como el número de repeticiones, se crea un control de entrenamiento mediante la función trainControl(), que se pasa al argumento trControl de la función train(). A efectos prácticos, cuando se aplican métodos de resampling hay que tener en cuenta dos cosas: el coste computacional que implica ajustar múltiples veces un modelo, cada vez con un subconjunto de datos distinto, y la reproducibilidad en la creación de las particiones. caret permite paralelizar el proceso para que sea más rápido y establecer semillas para asegurar que cada resampling o partición pueda crearse de nuevo con exactamente las mismas observaciones.

fitControl <- trainControl( method = "cv",  number = 10, savePredictions = 'final', classProbs = FALSE)

Lo siguiente será definir las variables que utilizaremos como predictoras y la variable a predecir:

predictors<-c("Sabor", "Color","Maltas")
outcomeName<-'Puntuacion'

Ahora solo nos falta correr el modelo utilizando la función train(). Lo que definirá qué modelo utilizaremos será el parámetro de method ya que la implementación es igual para todos los modelos en caret:

model1 <- train(train[,predictors], train[,outcomeName], method='lm', trControl = fitControl)
summary(model1)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -79.835  -2.084   1.231   3.274   8.635 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 83.59304    1.50626  55.497   <2e-16 ***
## Sabor       -0.02705    0.16635  -0.163    0.871    
## Color       -0.12969    0.15648  -0.829    0.408    
## Maltas       0.07947    0.11809   0.673    0.502    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.473 on 280 degrees of freedom
## Multiple R-squared:  0.003972,   Adjusted R-squared:  -0.006699 
## F-statistic: 0.3722 on 3 and 280 DF,  p-value: 0.7731

Una vez hemos obtenido los valores en el modelo entrenado podemos probar si los resultados son funcionales usando nuestro grupo de prueba. Para ello utilizaremos la función predict(), la cual nos permite para aplicar un modelo a un nuevo conjunto de datos

test$predicho <- predict(model1, test[,predictors])
plot(test$Puntuacion, test$predicho)

t.test(test$Puntuacion, test$predicho)
## 
##  Welch Two Sample t-test
## 
## data:  test$Puntuacion and test$predicho
## t = 0.29597, df = 67.929, p-value = 0.7682
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.9821891  1.3242820
## sample estimates:
## mean of x mean of y 
##  83.16706  82.99601

Ahora, lo interesante que tiene caret es que tiene mayor variedad, incluyendo backward and forward estimation, e incluso variaciones más robustas. En este link podrán encontrar todas las posibles variaciones: https://topepo.github.io/caret/train-models-by-tag.html#linear-regression

Regresiones logísticas

Ahora, ,miramemos un modelo más de regresión. la regresión logística es un tipo de análisis de regresión utilizado para predecir el resultado de una variable categórica en función de las variables independientes o predictoras. Es útil para modelar la probabilidad de un evento ocurriendo como función de otros factores.

La existencia de una relación significativa entre una variable cualitativa con dos niveles y una variable continua se puede estudiar mediante otros test estadísticos tales como t-test o ANOVA (un ANOVA de dos grupos es equivalente al t-test). Sin embargo, la regresión logística permite además calcular la probabilidad de que la variable dependiente pertenezca a cada una de las dos categorías en función del valor que adquiera la variable independiente. Supóngase que se quiere estudiar la relación entre los niveles de colesterol y los ataques de corazón. Para ello, se mide el colesterol de un grupo de personas y durante los siguientes 20 años se monitoriza que individuos han sufrido un ataque. Un t-test entre los niveles de colesterol de las personas que han sufrido ataque vs las que no lo han sufrido permitiría contrastar la hipótesis de que el colesterol y los ataques al corazón están asociados. Si además se desea conocer la probabilidad de que una persona con un determinado nivel de colesterol sufra un infarto en los próximos 20 años, o poder conocer cuánto tiene que reducir el colesterol un paciente para no superar un 50% de probabilidad de padecer un infarto en los próximos 20 años, se tiene que recurrir a la regresión logística.

La característica central de un modelo logístico es la siguiente: relaciona la probabilidad de un resultado con una función exponencial de una variable predictora. Al modelar la probabilidad de un resultado, un modelo logístico logra dos cosas. Primero, modela más directamente lo que nos interesa, que es una probabilidad o proporción, como la probabilidad de que un cliente determinado compre un producto o la proporción esperada de un segmento que responderá a una promoción. En segundo lugar, limita el modelo al rango apropiado para una proporción, que es [0, 1]. Un modelo lineal básico como se genera con lm() no tiene tal límite y podría estimar una probabilidad sin sentido como 1.05 o -0.04. La ecuación de la función logística es:

Figura 1: Función logística

Veamos un ejemplo de lo que hace la función aplicando la probabilidad logística con plogis():

plogis( exp (0) / (exp (0) + 1))
## [1] 0.6224593
plogis(-Inf)
## [1] 0
plogis ( -0.2)
## [1] 0.450166

Otra opción es utilizar la función sigmoidal, que se representa así:

\[\begin{equation} \text{función sigmoide} = \sigma(x) = \dfrac{1}{1 + e^{-x}} \tag{1} \end{equation}\]

Para valores de x muy grandes positivos, el valor de e^{-x} es aproximadamente 0 por lo que el valor de la función sigmoide es 1. Para valores de x muy grandes negativos, el valor e^{-x} tiende a infinito por lo que el valor de la función sigmoide es 0.

El análisis de regresión logística se enmarca en el conjunto de Modelos Lineales Generalizados (GLM por sus siglas en inglés). Las probabilidades que describen el posible resultado de un único ensayo se modelan, como una función de variables explicativas, utilizando una función logística. Por lo tanto, los modelos lineales generalizados se pueden utilizar para modelar recuentos de datos (como el número de compras) o intervalos de tiempo (como el tiempo que se pasa en un sitio web) o variables binarias (por ejemplo, compró / no compró). La característica común de todos los modelos GLM es que relacionan predictores distribuidos normalmente con un resultado no normal mediante una función conocida como enlace. Esto significa que pueden adaptarse a modelos para muchas distribuciones diferentes utilizando un marco único y coherente.

Regresión logística simple

La Regresión Logística Simple, desarrollada por David Cox en 1958, es un método de regresión que permite estimar la probabilidad de una variable cualitativa binaria en función de una variable cuantitativa. Una de las principales aplicaciones de la regresión logística es la de clasificación binaria, en el que las observaciones se clasifican en un grupo u otro dependiendo del valor que tome la variable empleada como predictor. Por ejemplo, clasificar a un individuo desconocido como hombre o mujer en función del tamaño de la mandíbula.

Es importante tener en cuenta que, aunque la regresión logística permite clasificar, se trata de un modelo de regresión que modela el logaritmo de la probabilidad de pertenecer a cada grupo. La asignación final se hace en función de las probabilidades predichas.

Para correr el modelo utilizaremos la misma sintaxis de la función lm() pero en este caso en su versión generalizada, como comentamos antes. Para lograrlo utilizaremos la función glm() utilizando como parámetro del modelo binomial:

reg1 <- glm(Compra ~ Puntuacion, data=data, family =  "binomial")

Una vez que hemos corrido nuestro modelo realizaremos la respectiva verificación utilizando la función summary():

summary(reg1)
## 
## Call:
## glm(formula = Compra ~ Puntuacion, family = "binomial", data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1723  -1.0773  -0.9369   1.2327   3.1819  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -5.23401    2.13904  -2.447   0.0144 *
## Puntuacion   0.05942    0.02559   2.322   0.0202 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 480.85  on 351  degrees of freedom
## Residual deviance: 474.26  on 350  degrees of freedom
## AIC: 478.26
## 
## Number of Fisher Scoring iterations: 4

¿Qué encontramos en este modelo? Primero, vemos la significancia de cada una de las variables. Para determinar la significancia individual de cada uno de los predictores introducidos en un modelo de regresión logística se emplea el estadístico Z y el test Wald chi-test.

Ahora, ¿qué significan los coeficientes de esta regresión? Para saber cómo funcionan estos. ¿Qué significa un coeficiente de 0.05? Podemos usarlo para calcular la asociación de los puntajes obtenidos con el factor de la compra del vino, examinando la proporción de éxito contra la proporción de no éxito. Para esto calcularemos los valores de la asociación con la función exp() y la función coef() sobre los coeficientes de la regresión:

exp(coef(reg1))
## (Intercept)  Puntuacion 
## 0.005332111 1.061223352

Esto muestra que el efecto de Puntuación es una razón de probabilidades estimada de 1.06, lo que significa que los jugadores tienen 1.06 veces más probabilidades de ser comprado cuando tienen una alta puntuación. Otra forma de pensar en esto es que la puntuación aumenta la probabilidad de jugar en un 6,1%.

A diferencia de los modelos de regresión lineal, en los modelos logísticos no existe un equivalente a R^2 que determine exactamente la varianza explicada por el modelo. Se han desarrollado diferentes métodos conocidos como pseudoR^2 que intentan aproximarse al concepto de R^2 pero que, aunque su rango oscila entre 0 y 1, no se pueden considerar equivalentes.

A la hora de evaluar la validez y calidad de un modelo de regresión logística, se analiza tanto el modelo en su conjunto como los predictores que lo forman. Se considera que el modelo es útil si es capaz de mostrar una mejora explicando las observaciones respecto al modelo nulo (sin predictores). El test Likelihood ratio calcula la significancia de la diferencia de residuos entre el modelo de interés y el modelo nulo. El estadístico sigue una distribución chi-cuadrado con grados de libertad equivalentes a la diferencia de grados de libertad de los dos modelos. Para ello utilizaremos la función anova():

anova(reg1, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Compra
## 
## Terms added sequentially (first to last)
## 
## 
##            Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
## NULL                         351     480.85           
## Puntuacion  1   6.5863       350     474.26  0.01028 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En este caso vemos que, comparado con el modelo nulo, la regresión es significativa.

Regresión logística múltiple

La regresión logística múltiple es una extensión de la regresión logística simple. Se basa en los mismos principios que la regresión logística simple (explicados anteriormente) pero ampliando el número de predictores. Los predictores pueden ser tanto continuos como categóricos (similar a la regresión lineal múltiple):

reg2 <- glm(Compra ~ Puntuacion + Maltas, data=data, family = "binomial")
summary(reg2)
## 
## Call:
## glm(formula = Compra ~ Puntuacion + Maltas, family = "binomial", 
##     data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6654  -0.3113  -0.1637   0.3657   3.3551  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.08581    3.42471  -1.485   0.1375    
## Puntuacion   0.10184    0.04178   2.437   0.0148 *  
## Maltas      -0.84440    0.07837 -10.774   <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: 480.85  on 351  degrees of freedom
## Residual deviance: 190.19  on 349  degrees of freedom
## AIC: 196.19
## 
## Number of Fisher Scoring iterations: 6

Podemos hacer las mismas revisiones sobre los betas para entender su comportamiento:

exp(coef(reg2))
## (Intercept)  Puntuacion      Maltas 
## 0.006183847 1.107202022 0.429815863

Además, podemos evaluar el modelo usando anova():

anova(reg2, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Compra
## 
## Terms added sequentially (first to last)
## 
## 
##            Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                         351     480.85             
## Puntuacion  1    6.586       350     474.26  0.01028 *  
## Maltas      1  284.073       349     190.19  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Otra forma de obtener esta “significancia” es a través del test LR. El test Likelihood ratio calcula la significancia de la diferencia de residuos entre el modelo de interés y el modelo nulo. El estadístico sigue una distribución chi-cuadrado con grados de libertad equivalentes a la diferencia de grados de libertad de los dos modelos. La forma de calcularlo es la siguiente:

dif_residuos <- reg2$null.deviance - reg2$deviance
df <- reg2$df.null - reg2$df.residual
p_value <- pchisq(q = dif_residuos,df = df, lower.tail = FALSE)
p_value
## [1] 7.657451e-64

En este caso el modelo en su conjunto es significativo.

Convertir probabilidad en clasificación

Una de las principales aplicaciones de un modelo de regresión logística es clasificar la variable cualitativa en función de valor que tome el predictor. Para conseguir esta clasificación, es necesario establecer un threshold de probabilidad a partir de la cual se considera que la variable pertenece a uno de los niveles.Para ello haremos una conversión en nuestros datos predichos usando la función ifelse():

predicciones <- ifelse(test = reg2$fitted.values > 0.5, yes = 1, no = 0)
t1 <- table(data$Compra, predicciones)
t1
##    predicciones
##       0   1
##   0 182  19
##   1  19 132

Esta tabla nos da información muy general sobre las variables a predecir, y se denomina matriz de confusión.

Figura 2: Matriz de confusión

Como podrán apreciar en la tabla, cada una de las posiciones que se obtienen en la matriz de confusión tiene un significado. Lo que nos gustaría saber son ciertas medidas que podemos calcular con la matriz. Estos serán nuestros indicadores de éxito en nuestra regresión logística:

Figura 3: Medidas de desempeño

Estas 4 medidas son solo algunas de aquellas que podremos calcular. Sin embargo, este análisis lo llevaremos a cabo cuando realicemos modelos de clasificación. Revisemos algunas de ellas:

182/(182+19)
## [1] 0.9054726
132/(132+19)
## [1] 0.8741722
(182+132)/352
## [1] 0.8920455

Esto nos muestra que nuestra regresión es bastante robusta y la podemos utilizar para predecir el resultado de una variable categórica binaria

Paquete CARET

Si queremos aplicar el paquete caret solamente haremos unas pequeñas variaciones. Primero hacemos nuestra partición:

trainIndex1 <- createDataPartition(data$Compra, p=0.8, list=F)
train1 <- data[trainIndex1,]
test1 <- data[-trainIndex1,]

Para ahorrar algo de tiempo, en esta ocasión usaremos los mismos parámetros de control que vimos al principio, por lo que ahora escogeremos las variables predictoras y predicha:

predictors1<-c("Puntuacion","Maltas")
outcomeName1<-'Compra'

Ahora deberemos correr el modelo modificando el método para que sea una regresión logística:

train1$Compra <- as.factor(train1$Compra)
model2 <- train(train1[,predictors1], train1[,outcomeName1], method='glm', trControl = fitControl)
summary(model2)
## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6770  -0.3616  -0.2022   0.3758   3.3434  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -5.10305    3.93875  -1.296   0.1951    
## Puntuacion   0.10068    0.04788   2.103   0.0355 *  
## Maltas      -0.78432    0.07990  -9.816   <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: 385.24  on 281  degrees of freedom
## Residual deviance: 163.20  on 279  degrees of freedom
## AIC: 169.2
## 
## Number of Fisher Scoring iterations: 6

Una vez hemos obtenido los valores en el modelo entrenado podemos probar si los resultados son funcionales usando nuestro grupo de prueba. Una vez más utilizaremos la función predict(), la cual nos permite para aplicar un modelo a un nuevo conjunto de datos

test1$Compra <- as.factor(test1$Compra)
test1$predicho <- predict(model2, test1[,predictors1])
plot(test1$Compra, test1$predicho)

confusionMatrix(test1$Compra, test1$predicho)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 34  6
##          1  1 29
##                                           
##                Accuracy : 0.9             
##                  95% CI : (0.8048, 0.9588)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : 1.138e-12       
##                                           
##                   Kappa : 0.8             
##                                           
##  Mcnemar's Test P-Value : 0.1306          
##                                           
##             Sensitivity : 0.9714          
##             Specificity : 0.8286          
##          Pos Pred Value : 0.8500          
##          Neg Pred Value : 0.9667          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4857          
##    Detection Prevalence : 0.5714          
##       Balanced Accuracy : 0.9000          
##                                           
##        'Positive' Class : 0               
##

Regresiones polinómicas

El análisis de regresión implica identificar relaciones entre variables. Dado un conjunto de variables independientes o predictoras, nuestro objetivo es predecir una variable dependiente o objetivo. En el análisis de regresión, la variable objetivo es numérica y las variables predictoras pueden ser numéricas, categóricas u ordinales.

Como vimos en las sesiones anteriores tenemos a nuestra disposición una serie de regresiones que pueden realizar amplias estimaciones sobre el comportamiento de nuestros dataframe, sin embargo en algunas ocasiones se aprecia que una regresión no retorna una aproximación satisfactoria. Algunas veces se debe a que los datos tienen un comportamiento difícil de predecir, sin embargo en otras ocasiones se aprecia que tienen un comportamiento diferente a los comunes pero con un patrón. En dichos casos se recomienda usar una regresión polinomica.

La regresión polinómica es un caso especial de regresión lineal en la que la relación entre los datos x e y son modelados usando un polinomio en vez de una linea. Puede ser usada cuando la relación entre las variables es no lineal. En el caso de R es considerada un caso especial de regresión linear múltiple al efectuarse con el comando lm().

Para hacer una regresión polinómica se usa la función lm , sin embargo se debe tener una consideración adicional. Como el modelo que buscamos solo requiere de una variable, se deben realizar ajustes sobre los parámetros que nos interesan, dichos ajustes pueden tener dos enfoques, los cuales analizaremos acá.

Ajustes basados en la creación manual de términos

Si intentamos usar en R nuestro parámetro de interés al cuadrado, cubo o una potencia diferente, haciendo la operación DENTRO DE LA FUNCIÓN, R no efectuará el cambio, por lo que estaremos suministrándole la misma variable más de una vez.

En este caso, si queremos establecer directamente el grado y los términos que deben aparecer en nuestra regresión, podemos crear un grupo de variables auxiliares que sean la variable de interés al cuadrado, cubo o potencia que se esté trabajando.

var1 <- (data$Votos)^2
reg3 <- lm(Puntuacion ~ var1, data = data)
summary(reg3)
## 
## Call:
## lm(formula = Puntuacion ~ var1, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -79.825  -2.338   1.139   3.066  10.415 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.282e+01  3.333e-01 248.530  < 2e-16 ***
## var1        2.123e-06  8.019e-07   2.647  0.00848 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.101 on 350 degrees of freedom
## Multiple R-squared:  0.01963,    Adjusted R-squared:  0.01683 
## F-statistic: 7.009 on 1 and 350 DF,  p-value: 0.008479

Ahora comparemos con la versión dentro del código:

reg4 <- lm(Puntuacion ~ (Votos) ^ 2, data=data)
summary(reg4)
## 
## Call:
## lm(formula = Puntuacion ~ (Votos)^2, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -79.377  -2.050   0.878   3.053  10.843 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 82.371819   0.353513 233.010  < 2e-16 ***
## Votos        0.005047   0.001173   4.303 2.19e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.005 on 350 degrees of freedom
## Multiple R-squared:  0.05025,    Adjusted R-squared:  0.04753 
## F-statistic: 18.52 on 1 and 350 DF,  p-value: 2.189e-05

La creación de dichas variables auxiliares depende de que potencia se desee trabajar para el modelo. Una vez creadas las variables auxiliares se usa el comando lm con normalidad, siendo el primer parametro los valores y, anexando la variable original como segundo parámetro y las demás variables como parámetros adicionales haciendo uso del operador ‘+’.

reg3 <- lm(Puntuacion ~ Votos + var1, data = data)
summary(reg3)
## 
## Call:
## lm(formula = Puntuacion ~ Votos + var1, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -78.954  -1.873   0.581   2.961  11.238 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.194e+01  3.806e-01 215.301  < 2e-16 ***
## Votos        1.217e-02  2.736e-03   4.446 1.18e-05 ***
## var1        -5.289e-06  1.841e-06  -2.873  0.00431 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.944 on 349 degrees of freedom
## Multiple R-squared:  0.07219,    Adjusted R-squared:  0.06687 
## F-statistic: 13.58 on 2 and 349 DF,  p-value: 2.097e-06
var2 <- (data$Votos)^(1/2)
reg3 <- lm(Puntuacion ~ Votos + var2, data = data)
summary(reg3)
## 
## Call:
## lm(formula = Puntuacion ~ Votos + var2, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -77.293  -1.722   0.413   2.674  12.287 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 79.736766   0.630279 126.510  < 2e-16 ***
## Votos       -0.010314   0.003288  -3.137  0.00185 ** 
## var2         0.566875   0.113882   4.978 1.01e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.811 on 349 degrees of freedom
## Multiple R-squared:  0.1132, Adjusted R-squared:  0.1081 
## F-statistic: 22.28 on 2 and 349 DF,  p-value: 7.854e-10
var3 <- (data$Votos)^(3)
reg3 <- lm(Puntuacion ~ Votos + var3, data = data)
summary(reg3)
## 
## Call:
## lm(formula = Puntuacion ~ Votos + var3, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -79.122  -1.876   0.661   3.035  11.085 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.211e+01  3.684e-01 222.895  < 2e-16 ***
## Votos        8.481e-03  1.878e-03   4.517 8.61e-06 ***
## var3        -1.614e-09  6.920e-10  -2.332   0.0203 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.968 on 349 degrees of freedom
## Multiple R-squared:  0.06482,    Adjusted R-squared:  0.05946 
## F-statistic:  12.1 on 2 and 349 DF,  p-value: 8.34e-06

Ajustes basados en I() yPoly()

En este ajuste lo que hacemos es decirle a R mediante la función I() y poly() que queremos realizar una regresión polinomica.

Usamos I() para incluir un término específico en un modelo, por ejemplo, I(x^3) indica que queremos solo el término X^{3}:

reg3 <- lm(Puntuacion ~ Votos + I(Votos^3), data = data)
summary(reg3)
## 
## Call:
## lm(formula = Puntuacion ~ Votos + I(Votos^3), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -79.122  -1.876   0.661   3.035  11.085 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.211e+01  3.684e-01 222.895  < 2e-16 ***
## Votos        8.481e-03  1.878e-03   4.517 8.61e-06 ***
## I(Votos^3)  -1.614e-09  6.920e-10  -2.332   0.0203 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.968 on 349 degrees of freedom
## Multiple R-squared:  0.06482,    Adjusted R-squared:  0.05946 
## F-statistic:  12.1 on 2 and 349 DF,  p-value: 8.34e-06

La función poly recibe 3 parametros: el primero es la variable de interés, la segunda es el grado del polinomio a crear y la tercera es “raw”, si raw está asignada con un valor verdadero la regresión se llevará a cabo con polinomios ajustados, si se configura como falsa se hará uso de polinomios ortogonales, proceso que resulta útil en espacios de trabajo diferentes al nuestro.

Si deseamos usar poly() podemos hacerlo dentro de nuestro comando de regresión, se ingresa como segundo parámetro directamente:

reg4 <- lm(Puntuacion ~  poly(Votos, 2) , data=data)
summary(reg4)
## 
## Call:
## lm(formula = Puntuacion ~ poly(Votos, 2), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -78.954  -1.873   0.581   2.961  11.238 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      83.0175     0.3168 262.036  < 2e-16 ***
## poly(Votos, 2)1  25.8413     5.9440   4.347 1.81e-05 ***
## poly(Votos, 2)2 -17.0776     5.9440  -2.873  0.00431 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.944 on 349 degrees of freedom
## Multiple R-squared:  0.07219,    Adjusted R-squared:  0.06687 
## F-statistic: 13.58 on 2 and 349 DF,  p-value: 2.097e-06