R Language
Modelos Lineales (Regresión)

Sintaxis

• lm (fórmula, datos, subconjunto, pesos, na.action, method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrastes = NULL, offset, .. .)

Parámetros

Regresión lineal en el conjunto de datos mtcars

El marco de datos mtcars incorporado contiene información sobre 32 automóviles, incluido su peso, eficiencia de combustible (en millas por galón), velocidad, etc. (Para obtener más información sobre el conjunto de datos, use la help(mtcars) ).

Si estamos interesados ​​en la relación entre la eficiencia del combustible ( mpg ) y el peso ( wt ), podemos comenzar a trazar esas variables con:

plot(mpg ~ wt, data = mtcars, col=2)

Las gráficas muestran una relación (lineal) !. Luego, si queremos realizar una regresión lineal para determinar los coeficientes de un modelo lineal, lm función lm :

fit <- lm(mpg ~ wt, data = mtcars)

El ~ aquí significa "explicado por", por lo que la fórmula mpg ~ wt significa que estamos prediciendo mpg como se explica por wt. La forma más útil de ver la salida es con:

summary(fit)

Lo que da la salida:

Call:
lm(formula = mpg ~ wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5432 -2.3647 -0.1252  1.4096  6.8727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Esto proporciona información sobre:

• la pendiente estimada de cada coeficiente ( wt y la intersección en y), que sugiere que la mejor predicción de mpg es 37.2851 + (-5.3445) * wt
• El valor p de cada coeficiente, lo que sugiere que la intersección y el peso probablemente no se deben al azar.
• Estimaciones generales de ajuste, tales como R ^ 2 y R ^ 2 ajustado, que muestran cuánta variación de mpg se explica por el modelo

Podríamos agregar una línea a nuestro primer gráfico para mostrar el mpg previsto:

abline(fit,col=3,lwd=2)

También es posible agregar la ecuación a esa gráfica. En primer lugar, obtener los coeficientes con coef . Luego, utilizando paste0 , paste0 los coeficientes con las variables apropiadas y +/- , para construir la ecuación. Finalmente, lo agregamos a la trama usando mtext :

bs <- round(coef(fit), 3) 
lmlab <- paste0("mpg = ", bs[1],
             ifelse(sign(bs[2])==1, " + ", " - "), abs(bs[2]), " wt ")
mtext(lmlab, 3, line=-2) 

El resultado es:

Trazando La Regresión (base)

Continuando con el ejemplo de mtcars , aquí hay una manera simple de producir una gráfica de su regresión lineal que es potencialmente adecuada para la publicación.

Primero ajuste el modelo lineal y

fit <- lm(mpg ~ wt, data = mtcars)

Luego trace las dos variables de interés y agregue la línea de regresión dentro del dominio de definición:

plot(mtcars$wt,mtcars$mpg,pch=18, xlab = 'wt',ylab = 'mpg')
lines(c(min(mtcars$wt),max(mtcars$wt)),
as.numeric(predict(fit, data.frame(wt=c(min(mtcars$wt),max(mtcars$wt))))))

¡Casi allí! El último paso es agregar a la gráfica, la ecuación de regresión, el cuadrado y el coeficiente de correlación. Esto se hace usando la función de vector :

rp = vector('expression',3)
rp[1] = substitute(expression(italic(y) == MYOTHERVALUE3 + MYOTHERVALUE4 %*% x), 
          list(MYOTHERVALUE3 = format(fit$coefficients[1], digits = 2),
                        MYOTHERVALUE4 = format(fit$coefficients[2], digits = 2)))[2]
rp[2] = substitute(expression(italic(R)^2 == MYVALUE), 
             list(MYVALUE = format(summary(fit)$adj.r.squared,dig=3)))[2]
rp[3] = substitute(expression(Pearson-R == MYOTHERVALUE2), 
             list(MYOTHERVALUE2 = format(cor(mtcars$wt,mtcars$mpg), digits = 2)))[2]

legend("topright", legend = rp, bty = 'n')

Tenga en cuenta que puede agregar cualquier otro parámetro como el RMSE adaptando la función de vector. Imagina que quieres una leyenda con 10 elementos. La definición del vector sería la siguiente:

rp = vector('expression',10)

y necesitarás definir r[1] .... a r[10]

Aquí está la salida:

Ponderación

A veces queremos que el modelo le dé más peso a algunos puntos de datos o ejemplos que a otros. Esto es posible especificando el peso de los datos de entrada mientras se aprende el modelo. En general, hay dos tipos de escenarios en los que podríamos usar ponderaciones no uniformes en los ejemplos:

• Pesos analíticos: Reflejan los diferentes niveles de precisión de diferentes observaciones. Por ejemplo, si el análisis de datos donde cada observación es el promedio de los resultados de un área geográfica, el peso analítico es proporcional a la inversa de la varianza estimada. Es útil cuando se trata de promedios en los datos al proporcionar un peso proporcional dado el número de observaciones. Fuente

• Pesos de muestreo (Pesos de probabilidad inversa - IPW): una técnica estadística para calcular estadísticas estandarizadas para una población diferente de aquella en la que se recopilaron los datos. Los diseños de estudio con una población de muestreo dispar y una población de inferencia objetivo (población objetivo) son comunes en la aplicación. Útil cuando se trata de datos que tienen valores perdidos. Fuente

La función lm() hace ponderación analítica. Para las ponderaciones de muestreo, el paquete de survey se usa para construir un objeto de diseño de encuesta y ejecutar svyglm() . Por defecto, el paquete de la survey utiliza ponderaciones muestrales. (NOTA: lm() y svyglm() con la familia gaussian() producirán las mismas estimaciones puntuales, ya que ambas resuelven los coeficientes minimizando los mínimos cuadrados ponderados. Difieren en cómo se calculan los errores estándar.)

Datos de prueba

data <- structure(list(lexptot = c(9.1595012302023, 9.86330744180814, 
8.92372556833205, 8.58202430280175, 10.1133857229336), progvillm = c(1L, 
1L, 1L, 1L, 0L), sexhead = c(1L, 1L, 0L, 1L, 1L), agehead = c(79L, 
43L, 52L, 48L, 35L), weight = c(1.04273509979248, 1.01139605045319, 
1.01139605045319, 1.01139605045319, 0.76305216550827)), .Names = c("lexptot", 
"progvillm", "sexhead", "agehead", "weight"), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -5L))

Pesos analíticos

lm.analytic <- lm(lexptot ~ progvillm + sexhead + agehead, 
                            data = data, weight = weight)
summary(lm.analytic)

Salida

Call:
lm(formula = lexptot ~ progvillm + sexhead + agehead, data = data, 
    weights = weight)

Weighted Residuals:
         1          2          3          4          5 
 9.249e-02  5.823e-01  0.000e+00 -6.762e-01 -1.527e-16

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.016054   1.744293   5.742    0.110
progvillm   -0.781204   1.344974  -0.581    0.665
sexhead      0.306742   1.040625   0.295    0.818
agehead     -0.005983   0.032024  -0.187    0.882

Residual standard error: 0.8971 on 1 degrees of freedom
Multiple R-squared:  0.467, Adjusted R-squared:  -1.132 
F-statistic: 0.2921 on 3 and 1 DF,  p-value: 0.8386

Pesos de muestreo (IPW)

library(survey)
data$X <- 1:nrow(data)             # Create unique id
        
# Build survey design object with unique id, ipw, and data.frame
des1 <- svydesign(id = ~X,  weights = ~weight, data = data)
        
# Run glm with survey design object
prog.lm <- svyglm(lexptot ~ progvillm + sexhead + agehead, design=des1)

Salida

Call:
svyglm(formula = lexptot ~ progvillm + sexhead + agehead, design = des1)

Survey design:
svydesign(id = ~X, weights = ~weight, data = data)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 10.016054   0.183942  54.452   0.0117 *
progvillm   -0.781204   0.640372  -1.220   0.4371  
sexhead      0.306742   0.397089   0.772   0.5813  
agehead     -0.005983   0.014747  -0.406   0.7546  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.2078647)

Number of Fisher Scoring iterations: 2

Comprobación de no linealidad con regresión polinomial

A veces, cuando se trabaja con regresión lineal, debemos verificar la no linealidad en los datos. Una forma de hacer esto es ajustar un modelo polinomial y verificar si se ajusta a los datos mejor que un modelo lineal. Existen otras razones, como las teóricas, que indican que se ajusta a un modelo cuadrático o de orden superior porque se cree que la relación de las variables es de naturaleza inherentemente polinomial.

mtcars un modelo cuadrático para el conjunto de datos mtcars . Para un modelo lineal, consulte Regresión lineal en el conjunto de datos mtcars .

Primero hacemos un diagrama de dispersión de las variables mpg (Millas / galón), disp (Desplazamiento (cu.in.)) y wt (Peso (1000 lbs)). La relación entre mpg y disp parece no lineal.

plot(mtcars[,c("mpg","disp","wt")])

Un ajuste lineal mostrará que disp no es significativo.

fit0 = lm(mpg ~ wt+disp, mtcars)
summary(fit0)

# Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 34.96055    2.16454  16.151 4.91e-16 ***
#wt          -3.35082    1.16413  -2.878  0.00743 ** 
#disp        -0.01773    0.00919  -1.929  0.06362 .  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.917 on 29 degrees of freedom
#Multiple R-squared:  0.7809,    Adjusted R-squared:  0.7658

Luego, para obtener el resultado de un modelo cuadrático, agregamos I(disp^2) . El nuevo modelo parece mejor cuando se observa R^2 y todas las variables son significativas.

fit1 = lm(mpg ~ wt+disp+I(disp^2), mtcars)
summary(fit1)

# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 41.4019837  2.4266906  17.061  2.5e-16 ***
#wt          -3.4179165  0.9545642  -3.581 0.001278 ** 
#disp        -0.0823950  0.0182460  -4.516 0.000104 ***
#I(disp^2)    0.0001277  0.0000328   3.892 0.000561 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.391 on 28 degrees of freedom
#Multiple R-squared:  0.8578,    Adjusted R-squared:  0.8426 

Como tenemos tres variables, el modelo ajustado es una superficie representada por:

mpg = 41.4020-3.4179*wt-0.0824*disp+0.0001277*disp^2

Otra forma de especificar la regresión polinomial es usar poly con el parámetro raw=TRUE , de lo contrario, se considerarán los polinomios ortogonales (consulte la help(ploy) para obtener más información). Obtenemos el mismo resultado usando:

summary(lm(mpg ~ wt+poly(disp, 2, raw=TRUE),mtcars))

Por último, ¿qué pasa si necesitamos mostrar una gráfica de la superficie estimada? Bueno, hay muchas opciones para hacer gráficos 3D en R Aquí usamos Fit3d del paquete p3d .

library(p3d)
Init3d(family="serif", cex = 1)
Plot3d(mpg ~ disp+wt, mtcars)
Axes3d()
Fit3d(fit1)

Evaluación de la calidad

Después de construir un modelo de regresión, es importante verificar el resultado y decidir si el modelo es apropiado y funciona bien con los datos disponibles. Esto se puede hacer examinando el gráfico de residuos, así como otros gráficos de diagnóstico.

# fit the model
fit <- lm(mpg ~ wt, data = mtcars)
# 
par(mfrow=c(2,1))
# plot model object
plot(fit, which =1:2)

Estas parcelas verifican dos suposiciones que se hicieron al construir el modelo:

1. Que el valor esperado de la variable predicha (en este caso mpg ) viene dado por una combinación lineal de los predictores (en este caso wt ). Esperamos que esta estimación sea imparcial. Por lo tanto, los residuos deben centrarse alrededor de la media de todos los valores de los predictores. En este caso, vemos que los residuos tienden a ser positivos en los extremos y negativos en el medio, lo que sugiere una relación no lineal entre las variables.
2. Que la variable predicha real se distribuye normalmente alrededor de su estimación. Por lo tanto, los residuos deben ser distribuidos normalmente. Para los datos distribuidos normalmente, los puntos en una gráfica QQ normal deben estar en o cerca de la diagonal. Hay una cierta cantidad de sesgo en los extremos aquí.

Usando la función 'predecir'

Una vez que se construye un modelo, predict es la función principal para probar con nuevos datos. Nuestro ejemplo utilizará el conjunto de datos incorporado de mtcars para retroceder millas por galón contra desplazamiento:

my_mdl <- lm(mpg ~ disp, data=mtcars)
my_mdl

Call:
lm(formula = mpg ~ disp, data = mtcars)

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

Si tuviera una nueva fuente de datos con desplazamiento, podría ver las millas estimadas por galón.

set.seed(1234)
newdata <- sample(mtcars$disp, 5)
newdata
[1] 258.0  71.1  75.7 145.0 400.0

newdf <- data.frame(disp=newdata)
predict(my_mdl, newdf)
       1        2        3        4        5 
18.96635 26.66946 26.47987 23.62366 13.11381

La parte más importante del proceso es crear un nuevo marco de datos con los mismos nombres de columna que los datos originales. En este caso, los datos originales tenían una columna etiquetada disp , estaba seguro de llamar a los nuevos datos con el mismo nombre.

Precaución

Echemos un vistazo a algunas trampas comunes:

1. no usar un data.frame en el nuevo objeto:

   predict(my_mdl, newdata)
   Error in eval(predvars, data, env) : 
      numeric 'envir' arg not of length one
   

2. no usar los mismos nombres en el nuevo marco de datos:

   newdf2 <- data.frame(newdata)
   predict(my_mdl, newdf2)
   Error in eval(expr, envir, enclos) : object 'disp' not found
   

Exactitud

Para verificar la precisión de la predicción, necesitará los valores y reales de los nuevos datos. En este ejemplo, newdf necesitará una columna para 'mpg' y 'disp'.

newdf <- data.frame(mpg=mtcars$mpg[1:10], disp=mtcars$disp[1:10])
#     mpg  disp
# 1  21.0 160.0
# 2  21.0 160.0
# 3  22.8 108.0
# 4  21.4 258.0
# 5  18.7 360.0
# 6  18.1 225.0
# 7  14.3 360.0
# 8  24.4 146.7
# 9  22.8 140.8
# 10 19.2 167.6

p <- predict(my_mdl, newdf)

#root mean square error
sqrt(mean((p - newdf$mpg)^2, na.rm=TRUE))
[1] 2.325148