# TRABAJO FINAL - REGRESION LINEAL MULTIPLE
# NASA Airfoil Self-Noise
# Integrantes: Miguel Turizo, Carlos Garay, Angel Mendoza, Valentina Puentes
library(readxl)
library(nortest)
library(lmtest)
library(car)
# Si NASA.xlsx esta en otra carpeta, sustituir solo esta ruta.
NASA <- read_excel("NASA.xlsx", sheet = "Base de datos")
SSPL <- NASA$SSPL
f <- NASA$f
alpha <- NASA$alpha
c <- NASA$c
U_infinity <- NASA$U_infinity
delta <- NASA$delta
n <- length(SSPL)
k <- 5
Universidad Pontificia Bolivariana
Regresión lineal, semestre 2026-20
Modelo de regresión lineal múltiple del ruido aerodinámico
Integrantes
Carlos José Garay Hoyos
Ángel Antonio Mendoza Burgos
Valentina Puentes Vanegas
Miguel Andrés Turizo Barrios
Profesor
Delio Salgado Chamorro
Fecha de entrega
10 de noviembre de 2026
Introducción
Cuando el aire fluye sobre el ala de un avión, una pala de helicóptero, una turbina eólica o un ventilador, la propia superficie genera ruido. Este ruido autogenerado del perfil aerodinámico (airfoil self-noise) se produce sobre todo porque la capa de aire pegada a la superficie (la capa límite) se vuelve turbulenta e interactúa con el borde de salida del perfil. Reducirlo importa para el confort de las comunidades cercanas a aeropuertos y parques eólicos y para cumplir normas de ruido.
La base proviene de ensayos de la NASA en un túnel de viento anecoico (sin ecos), con perfiles de distintos tamaños, a distintas velocidades y ángulos de ataque. Cada fila es una medición del ruido a una frecuencia dada bajo unas condiciones de ensayo. Objetivo: construir un modelo de regresión que explique y prediga ese ruido a partir de las condiciones de operación y verificar si cumple los supuestos de la regresión lineal.
Este informe analiza el nivel de presión sonora SSPL registrado en la base NASA Airfoil Self-Noise. Se estudia su relación con frecuencia (f), ángulo de ataque (alpha), longitud de cuerda (c), velocidad del flujo (U_infinity) y espesor de desplazamiento (delta). El menú lateral y la tabla siguiente permiten ir directamente a cada apartado. Los bloques de código aparecen completos y, debajo de ellos, se explica cómo leer sus salidas.
Objetivo: estimar la relación entre las cinco características aerodinámicas y el nivel de presión sonora, comprobar los supuestos del modelo y comparar alternativas para seleccionar una especificación final.
Tabla de contenido
- Introducción
- Datos y preparación
- Exploración
- Modelo y estimación de parámetros
- Inferencia sobre el modelo completo
- Inferencia para cada parámetro
- Normalidad de los residuales
- Varianza constante
- Independencia
- Multicolinealidad
- Observaciones influyentes
- Selección: todas las regresiones posibles
- Métodos paso a paso
- Modelo elegido y conclusión
- Revisión de supuestos del modelo final
- Cierre
Datos y preparación
La información se importa desde la hoja Base de datos de NASA.xlsx. Guarda el Excel y este archivo R Markdown en la misma carpeta antes de usar Knit. Las asignaciones con $ extraen las columnas; n cuenta las observaciones y k fija las cinco regresoras del modelo original.
| Variable | Papel en el análisis | Unidad |
|---|---|---|
SSPL |
Nivel de presión sonora, respuesta | dB |
f |
Frecuencia | Hz |
alpha |
Ángulo de ataque | grados |
c |
Longitud de cuerda | m |
U_infinity |
Velocidad del flujo | m/s |
delta |
Espesor de desplazamiento | m |
La variable respuesta es SSPL (dB). f se mide en Hz, alpha en grados, c en metros, U_infinity en m/s y delta en metros. La base importada contiene 1503 filas y 0 valores faltantes.
Exploracion
summary(NASA) presenta mínimos, cuartiles, medianas, medias y máximos. colSums(is.na(NASA)) comprueba datos faltantes; cor(NASA) calcula correlaciones lineales y pairs(NASA) dibuja una matriz de dispersión.
## Exploracion
summary(NASA) f alpha c U_infinity
Min. : 200 Min. : 0.000 Min. :0.0254 Min. :31.70
1st Qu.: 800 1st Qu.: 2.000 1st Qu.:0.0508 1st Qu.:39.60
Median : 1600 Median : 5.400 Median :0.1016 Median :39.60
Mean : 2886 Mean : 6.782 Mean :0.1365 Mean :50.86
3rd Qu.: 4000 3rd Qu.: 9.900 3rd Qu.:0.2286 3rd Qu.:71.30
Max. :20000 Max. :22.200 Max. :0.3048 Max. :71.30
delta SSPL
Min. :0.0004007 Min. :103.4
1st Qu.:0.0025351 1st Qu.:120.2
Median :0.0049574 Median :125.7
Mean :0.0111399 Mean :124.8
3rd Qu.:0.0155759 3rd Qu.:130.0
Max. :0.0584113 Max. :141.0 colSums(is.na(NASA)) f alpha c U_infinity delta SSPL
0 0 0 0 0 0 round(cor(NASA), 3) f alpha c U_infinity delta SSPL
f 1.000 -0.273 -0.004 0.134 -0.230 -0.391
alpha -0.273 1.000 -0.505 0.059 0.753 -0.156
c -0.004 -0.505 1.000 0.004 -0.221 -0.236
U_infinity 0.134 0.059 0.004 1.000 -0.004 0.125
delta -0.230 0.753 -0.221 -0.004 1.000 -0.313
SSPL -0.391 -0.156 -0.236 0.125 -0.313 1.000pairs(NASA, pch = 20, cex = 0.4)
Análisis de la gráfica: cada panel compara dos variables. Una nube ascendente sugiere relación positiva; una descendente, negativa; una nube curva indica que una recta podría no describir bien esa pareja. Las correlaciones y estos gráficos son exploratorios: no equivalen a los efectos parciales del modelo múltiple.
1. Modelo y estimacion de parametros
lm() ajusta por mínimos cuadrados el modelo completo. summary(modelo) muestra los coeficientes estimados, errores estándar, estadísticos t, valores p y medidas de ajuste.
Fórmulas: 1. Modelo y estimación de parámetros
Modelo poblacional (con \(k = 5\) regresoras):
\[ \text{SSPL}_i = \beta_0 + \beta_1 f_i + \beta_2\,\text{alpha}_i + \beta_3 c_i + \beta_4 U_{\infty,i} + \beta_5\,\delta_i + \varepsilon_i, \qquad \varepsilon_i \sim N(0,\sigma^2)\ \text{independientes} \]
Forma matricial:
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]
Estimador de mínimos cuadrados:
\[ \hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y} \]
Valores ajustados, residuales y suma de cuadrados del error:
\[ \hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\beta}}, \qquad \mathbf{e} = \mathbf{y}-\hat{\mathbf{y}}, \qquad SSE = \sum_{i=1}^{n} e_i^2 \]
Estimador de la varianza y errores estándar:
\[ \hat{\sigma}^2 = MSE = \frac{SSE}{n-k-1}, \qquad \widehat{Var}(\hat{\boldsymbol{\beta}}) = MSE\,(\mathbf{X}'\mathbf{X})^{-1}, \qquad EE(\hat{\beta}_j)=\sqrt{MSE\,c_{jj}} \]
donde \(c_{jj}\) es el \(j\)-ésimo elemento de la diagonal de \((\mathbf{X}'\mathbf{X})^{-1}\).
Ecuación estimada del modelo original:
\[ \widehat{\text{SSPL}} = 132.8 - 0.001282\,f - 0.4219\,\text{alpha} - 35.69\,c + 0.09985\,U_\infty - 147.3\,\delta \]
Coeficientes de determinación:
\[ R^2 = 1-\frac{SSE}{SST}, \qquad R^2_{aj} = 1-\frac{n-1}{n-k-1}\cdot\frac{SSE}{SST} \]
## 1. Modelo y estimacion de parametros
# SSPL = b0 + b1*f + b2*alpha + b3*c + b4*U_infinity + b5*delta + e
modelo <- lm(SSPL ~ f + alpha + c + U_infinity + delta, data = NASA)
summary(modelo)
Call:
lm(formula = SSPL ~ f + alpha + c + U_infinity + delta, data = NASA)
Residuals:
Min 1Q Median 3Q Max
-17.480 -2.882 -0.209 3.152 16.064
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.328e+02 5.447e-01 243.87 <2e-16 ***
f -1.282e-03 4.211e-05 -30.45 <2e-16 ***
alpha -4.219e-01 3.890e-02 -10.85 <2e-16 ***
c -3.569e+01 1.630e+00 -21.89 <2e-16 ***
U_infinity 9.985e-02 8.132e-03 12.28 <2e-16 ***
delta -1.473e+02 1.501e+01 -9.81 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.809 on 1497 degrees of freedom
Multiple R-squared: 0.5157, Adjusted R-squared: 0.5141
F-statistic: 318.8 on 5 and 1497 DF, p-value: < 2.2e-16La ecuación estimada combina un intercepto y una pendiente por regresora. Cada pendiente representa el cambio esperado en SSPL al aumentar una unidad esa variable manteniendo constantes las otras cuatro. El R² del modelo original es 0.5157 y el R² ajustado es 0.5141; estas medidas describen ajuste, no validez automática de los supuestos.
2. Inferencia sobre el modelo completo
modelo_nulo incluye solo el intercepto. anova(modelo_nulo, modelo) compara su error con el del modelo de cinco regresoras. La prueba F evalúa si, en conjunto, las variables aportan información sobre SSPL.
Fórmulas: 2. Inferencia sobre el modelo completo
Hipótesis:
\[ H_0:\ \beta_1=\beta_2=\beta_3=\beta_4=\beta_5=0 \qquad \text{vs} \qquad H_1:\ \text{al menos un } \beta_j \neq 0 \]
Descomposición de la variabilidad:
\[ SST = SSR + SSE, \qquad SST=\sum_{i=1}^{n}(y_i-\bar{y})^2, \qquad SSR=\sum_{i=1}^{n}(\hat{y}_i-\bar{y})^2 \]
Grados de libertad y cuadrados medios:
\[ gl_{SSR}=k,\quad gl_{SSE}=n-k-1,\quad gl_{SST}=n-1, \qquad MSR=\frac{SSR}{k}, \quad MSE=\frac{SSE}{n-k-1} \]
Estadístico de prueba y regla de decisión:
\[ F_0 = \frac{MSR}{MSE} \sim F_{k,\,n-k-1} \]
\[ \text{Se rechaza } H_0 \text{ si } F_0 > F_{\alpha,\,k,\,n-k-1} \quad \text{o si } \text{valor-}p = P(F_{k,n-k-1} > F_0) < \alpha \]
Prueba entre modelos anidados (modelo reducido R frente a completo C, con \(q\) coeficientes igualados a cero):
\[ F_0=\frac{(SSE_R-SSE_C)/q}{MSE_C}\sim F_{q,\,n-k-1} \]
## 2. Inferencia sobre el modelo completo
# H0: b1 = b2 = b3 = b4 = b5 = 0; H1: al menos uno es distinto de 0.
modelo_nulo <- lm(SSPL ~ 1, data = NASA)
anova(modelo_nulo, modelo)Analysis of Variance Table
Model 1: SSPL ~ 1
Model 2: SSPL ~ f + alpha + c + U_infinity + delta
Res.Df RSS Df Sum of Sq F Pr(>F)
1 1502 71482
2 1497 34618 5 36864 318.82 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1qf(0.05, k, n - k - 1, lower.tail = FALSE)[1] 2.220076El F observado es 318.82 y el valor crítico al 5 % es 2.22. La decisión es rechazar H₀: al menos una pendiente es distinta de cero. Esto no significa que todas las variables sean significativas por separado.
3. Inferencia para cada parametro
La tabla de summary(modelo)$coefficients permite estudiar un parámetro a la vez. Para cada pendiente, se contrasta cero frente a un valor distinto de cero; se compara |t| con el t crítico o se revisa si su valor p es menor que 0.05. confint() ofrece los intervalos al 95 %.
Fórmulas: 3. Inferencia para cada parámetro
Hipótesis (para \(j=0,1,\dots,5\)):
\[ H_0:\ \beta_j = 0 \qquad \text{vs} \qquad H_1:\ \beta_j \neq 0 \]
Error estándar y estadístico de prueba:
\[ EE(\hat{\beta}_j)=\sqrt{MSE\,[(\mathbf{X}'\mathbf{X})^{-1}]_{jj}}, \qquad t_0=\frac{\hat{\beta}_j}{EE(\hat{\beta}_j)} \sim t_{n-k-1} \]
Regla de decisión: se rechaza \(H_0\) si el estadístico de prueba supera al estadístico teórico, es decir, si
\[ |t_0| > t_{\alpha/2;\,n-k-1} \]
Intervalo de confianza del \((1-\alpha)\,100\,\%\):
\[ \hat{\beta}_j \pm t_{\alpha/2,\,n-k-1}\; EE(\hat{\beta}_j) \]
## 3. Inferencia para cada parametro
# H0: bj = 0; H1: bj es distinto de 0.
summary(modelo)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.328338e+02 5.447007e-01 243.865682 0.000000e+00
f -1.282207e-03 4.210547e-05 -30.452267 6.421817e-159
alpha -4.219117e-01 3.889609e-02 -10.847149 1.923376e-26
c -3.568800e+01 1.630432e+00 -21.888679 2.286839e-92
U_infinity 9.985404e-02 8.132259e-03 12.278758 4.340636e-33
delta -1.473005e+02 1.501467e+01 -9.810441 4.617244e-22t_teorico <- qt(0.025, df = n - k - 1, lower.tail = FALSE)
t_teorico[1] 1.96155confint(modelo, level = 0.95) 2.5 % 97.5 %
(Intercept) 1.317653e+02 1.339023e+02
f -1.364799e-03 -1.199615e-03
alpha -4.982083e-01 -3.456151e-01
c -3.888617e+01 -3.248983e+01
U_infinity 8.390221e-02 1.158059e-01
delta -1.767525e+02 -1.178485e+02Lectura: una pendiente con valor p menor que 0.05 presenta evidencia estadística de asociación condicional a las demás variables incluidas. El signo del coeficiente indica la dirección estimada; su magnitud debe leerse en las unidades de la regresora. No se debe confundir significancia estadística con importancia práctica.
4. Normalidad de los residuales
Se revisa si la distribución de los residuales es compatible con normalidad. El código guarda los residuales, traza un gráfico Q-Q, dibuja un histograma y ejecuta Lilliefors, Anderson-Darling y Shapiro-Wilk.
Fórmulas: 4. Normalidad de los residuales
Hipótesis (válidas para Lilliefors, Anderson-Darling y Shapiro-Wilk):
\[ H_0:\ \varepsilon_i \sim N(\mu,\sigma^2) \qquad \text{vs} \qquad H_1:\ \varepsilon_i \text{ no sigue una distribución normal} \]
Se rechaza \(H_0\) si el valor-\(p<\alpha\).
Con los residuales ordenados \(e_{(1)}\le\cdots\le e_{(n)}\) y \(z_{(i)}=(e_{(i)}-\bar{e})/s_e\), las tres pruebas usan los siguientes estadísticos.
Lilliefors (Kolmogorov-Smirnov con media y varianza estimadas):
\[ D=\max_i\left\{\frac{i}{n}-\Phi\!\left(z_{(i)}\right),\ \Phi\!\left(z_{(i)}\right)-\frac{i-1}{n}\right\} \]
Con \(n>30\) y \(\alpha=0.05\), el estadístico teórico aproximado es \(D_{teo}=0.886/\sqrt{n}\) (Lilliefors, 1967), y se rechaza \(H_0\) si el estadístico de prueba supera al teórico, es decir, si \(D>D_{teo}\).
Anderson-Darling:
\[ A^2=-n-\frac{1}{n}\sum_{i=1}^{n}(2i-1)\left[\ln\Phi(z_{(i)})+\ln\!\big(1-\Phi(z_{(n+1-i)})\big)\right], \qquad z_{(i)}=\frac{e_{(i)}-\bar{e}}{s} \]
Shapiro-Wilk:
\[ W=\frac{\left(\sum_{i=1}^{n}a_i\,e_{(i)}\right)^2}{\sum_{i=1}^{n}(e_i-\bar{e})^2} \]
Para Anderson-Darling y Shapiro-Wilk no hay un único estadístico teórico de tabla, por lo que la decisión se toma con el p-valor: se rechaza \(H_0\) si es menor que \(\alpha\).
## 4. Normalidad de los residuales
residuales <- residuals(modelo)
qqnorm(residuales, main = "Q-Q de residuales")
qqline(residuales, col = "red")
hist(residuales, breaks = 30, col = "steelblue",
main = "Histograma de residuales", xlab = "Residuales")
lillie.test(residuales)
Lilliefors (Kolmogorov-Smirnov) normality test
data: residuales
D = 0.027418, p-value = 0.01026ad.test(residuales)
Anderson-Darling normality test
data: residuales
A = 1.3911, p-value = 0.001339shapiro.test(residuales)
Shapiro-Wilk normality test
data: residuales
W = 0.99563, p-value = 0.0002465
El estadístico de prueba de Lilliefors y su estadístico teórico se calculan de la siguiente manera:
ord <- sort(residuales)
z <- (ord - mean(residuales)) / sd(residuales)
Fz <- pnorm(z)
i <- 1:n
D <- max(pmax(i / n - Fz, Fz - (i - 1) / n))
Dcrit <- 0.886 / sqrt(n)
data.frame(D_manual = D, D_teorico = Dcrit) D_manual D_teorico
1 0.0274179 0.02285358Para Lilliefors, el estadístico de prueba es D = 0.0274 y el estadístico teórico es 0.0229. Como D es mayor que el estadístico teórico, se rechaza H₀.
Análisis de la gráfica: en el Q-Q, puntos cercanos a la línea roja favorecen la aproximación normal; desviaciones sistemáticas, especialmente en las colas, la cuestionan. El histograma permite ver asimetría o colas largas. En las pruebas, un valor p inferior a 0.05 indica evidencia contra la normalidad; con 1503 datos, también pueden detectarse desviaciones pequeñas que deben contrastarse con los gráficos.
5. Varianza constante
fitted(modelo) obtiene los valores predichos. El diagrama de residuales frente a predichos y las dos versiones de Breusch-Pagan examinan si la dispersión del error cambia con el nivel estimado de SSPL.
Fórmulas: 5. Varianza constante
Hipótesis:
\[ H_0:\ Var(\varepsilon_i)=\sigma^2 \ \text{para todo } i \quad (\text{homocedasticidad}) \qquad \text{vs} \qquad H_1:\ Var(\varepsilon_i) \text{ depende de las regresoras} \]
Breusch-Pagan original (Cook-Weisberg, studentize = FALSE):
\[ \tilde{\sigma}^2=\frac{SSE}{n}, \qquad g_i=\frac{e_i^2}{\tilde{\sigma}^2}, \qquad BP=\tfrac{1}{2}\,SSR_{aux}\sim \chi^2_{k} \]
donde \(SSR_{aux}\) es la suma de cuadrados explicada al regresar \(g_i\) sobre las regresoras.
Breusch-Pagan studentizado (Koenker):
\[ BP_{s}=n\,R^2_{aux}\sim \chi^2_{k} \]
donde \(R^2_{aux}\) es el \(R^2\) al regresar \(e_i^2\) sobre las regresoras.
Regla de decisión:
\[ \text{Se rechaza } H_0 \text{ si } BP > \chi^2_{\alpha,\,k}\ \text{ o valor-}p<\alpha \]
## 5. Varianza constante
ajustados <- fitted(modelo)
plot(ajustados, residuales, pch = 20,
main = "Residuales vs ajustados",
xlab = "Valores ajustados", ylab = "Residuales")
abline(h = 0, col = "red")
bptest(modelo)
studentized Breusch-Pagan test
data: modelo
BP = 91.141, df = 5, p-value < 2.2e-16bptest(modelo, studentize = FALSE)
Breusch-Pagan test
data: modelo
BP = 116.32, df = 5, p-value < 2.2e-16qchisq(0.05, k, lower.tail = FALSE)[1] 11.0705Análisis de la gráfica: una banda de anchura aproximadamente similar alrededor de cero es compatible con varianza constante. Un embudo, una curva o zonas con dispersión creciente sugieren problemas. En Breusch-Pagan, valor p menor que 0.05 apoya heterocedasticidad. Las dos pruebas impresas usan versiones diferentes; conviene identificar cuál se cita al redactar la conclusión.
6. Independencia
El gráfico coloca los residuales en el orden de las filas para buscar rachas o patrones. durbinWatsonTest() calcula un estadístico para autocorrelación de primer orden y su prueba bilateral.
Fórmulas: 6. Independencia
Hipótesis (prueba bilateral), con \(\varepsilon_t=\rho\,\varepsilon_{t-1}+u_t\):
\[ H_0:\ \rho=0 \qquad \text{vs} \qquad H_1:\ \rho\neq 0 \]
Estadístico de Durbin-Watson:
\[ d=\frac{\sum_{t=2}^{n}(e_t-e_{t-1})^2}{\sum_{t=1}^{n}e_t^2}\approx 2\,(1-\hat{\rho}) \]
Regla de decisión (con las cotas \(d_L\) y \(d_U\) de la tabla):
| Valor de \(d\) | Decisión |
|---|---|
| \(d<d_L\) | Se rechaza \(H_0\): autocorrelación positiva |
| \(d_L\le d\le d_U\) | Prueba no concluyente |
| \(d_U<d<4-d_U\) | No se rechaza \(H_0\) |
| \(4-d_U\le d\le 4-d_L\) | Prueba no concluyente |
| \(d>4-d_L\) | Se rechaza \(H_0\): autocorrelación negativa |
Las tablas de Durbin-Watson no cubren \(n=1503\), así que el estadístico teórico también puede obtenerse con la aproximación normal: \(d\) tiene media 2 y varianza aproximada \(4/n\), por lo que se rechaza \(H_0\) si \(d\) queda fuera del intervalo
\[ 2 \pm z_{\alpha/2}\,\frac{2}{\sqrt{n}} \]
## 6. Independencia
plot(seq_len(n), residuales, pch = 20,
main = "Residuales segun el orden",
xlab = "Observacion", ylab = "Residuales")
abline(h = 0, col = "red")
set.seed(123)
durbinWatsonTest(modelo, alternative = "two.sided") lag Autocorrelation D-W Statistic p-value
1 0.775644 0.4474251 0
Alternative hypothesis: rho != 0# El orden de filas corresponde al experimento; no implica una serie temporal.El estadístico y los límites de no rechazo se calculan así:
d <- sum(diff(residuales)^2) / sum(residuales^2)
r1 <- sum(residuales[-1] * residuales[-n]) / sum(residuales^2)
z_crit <- qnorm(0.975)
d_inf <- 2 - z_crit * 2 / sqrt(n)
d_sup <- 2 + z_crit * 2 / sqrt(n)
data.frame(d = d, r1 = r1, limite_inferior = d_inf, limite_superior = d_sup) d r1 limite_inferior limite_superior
1 0.4474251 0.775644 1.898889 2.101111El estadístico de prueba es d = 0.4474, con una autocorrelación estimada de 0.776, y la región de no rechazo es el intervalo (1.8989; 2.1011). Como d cae fuera del intervalo, se rechaza H₀.
Análisis de la gráfica: alternancia sin patrón claro alrededor de cero favorece la independencia; tramos prolongados con el mismo signo sugieren dependencia o estructura omitida. Un Durbin-Watson cercano a 2 es compatible con ausencia de autocorrelación de primer orden. Cautela: los registros están ordenados por condiciones de ensayo, no por tiempo; el resultado no debe interpretarse como evidencia de una serie temporal.
7. Multicolinealidad
vif(modelo) informa cuánto aumenta la varianza estimada de cada pendiente debido a su relación con las otras regresoras.
Fórmulas: 7. Multicolinealidad
No es una prueba de hipótesis, sino un criterio de diagnóstico:
\[ VIF_j=\frac{1}{1-R_j^2} \]
donde \(R_j^2\) es el coeficiente de determinación de la regresión de \(x_j\) sobre las demás regresoras. Como criterio, un \(VIF<5\) indica que no hay problema, un \(VIF\) entre 5 y 10 indica multicolinealidad moderada y un \(VIF>10\) indica multicolinealidad severa.
## 7. Multicolinealidad
vif(modelo) f alpha c U_infinity delta
1.144444 3.441658 1.510754 1.041698 2.532127 VIF cercano a 1 indica poca redundancia lineal; valores elevados, por ejemplo superiores a 5 o 10 según el criterio empleado, ameritan revisar las relaciones entre regresoras. Un VIF alto puede volver inestables las pendientes individuales aunque la prueba global resulte significativa.
7.1. Observaciones influyentes
La distancia de Cook combina tamaño del residual e influencia sobre el ajuste. La línea roja marca 4/n como umbral exploratorio; el código también cuenta residuales estandarizados con valor absoluto mayor que 3.
Fórmulas: 7.1. Observaciones influyentes
\[ h_{ii}=\mathbf{x}_i'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i, \qquad r_i=\frac{e_i}{\sqrt{MSE\,(1-h_{ii})}}, \qquad D_i=\frac{r_i^2}{k+1}\cdot\frac{h_{ii}}{1-h_{ii}} \]
\[ \text{Observación influyente si } D_i>\frac{4}{n}, \qquad \text{observación atípica si } |r_i|>3 \]
## 7.1. Observaciones influyentes
cook <- cooks.distance(modelo)
plot(cook, type = "h", main = "Distancia de Cook",
xlab = "Observacion", ylab = "Distancia de Cook")
abline(h = 4 / n, col = "red")
sum(cook > 4 / n)[1] 114sum(abs(rstandard(modelo)) > 3)[1] 13Análisis de la gráfica: las barras que sobresalen del umbral son casos para inspeccionar en la base y contrastar con el experimento. No se eliminan automáticamente: un punto inusual puede ser una medición válida.
8. Seleccion del modelo: todas las regresiones posibles
Se forman logf_c = log(f) - media(log(f)), logf_c2 = logf_c² y logdelta = log(delta). La escala logarítmica comprime los valores grandes de f y delta; el término cuadrático permite que la relación con frecuencia sea curva. Centrar log(f) facilita interpretar el intercepto y reduce parte de la redundancia con su cuadrado. Son transformaciones de variables medidas, no datos nuevos.
Fórmulas: Mejora del modelo (transformaciones)
\[ \log f_c=\ln f-\frac{1}{n}\sum_{i=1}^{n}\ln f_i, \qquad \log f_{c}^{2}=(\log f_c)^2, \qquad \log\delta=\ln\delta \]
Modelo con las seis variables candidatas (\(k=6\)):
\[ \text{SSPL}=\beta_0+\beta_1\log f_c+\beta_2(\log f_c)^2+\beta_3\,\text{alpha}+\beta_4\,c+\beta_5\,U_\infty+\beta_6\log\delta+\varepsilon \]
Sus hipótesis son las mismas de los apartados 2 y 3, con \(k=6\):
\[ H_0:\ \beta_1=\dots=\beta_6=0 \quad\text{vs}\quad H_1:\ \text{algún }\beta_j\neq0; \qquad H_0:\ \beta_j=0 \quad\text{vs}\quad H_1:\ \beta_j\neq0 \]
El bucle recorre los 63 subconjuntos no vacíos de seis términos candidatos. Después se descartan los que incluyen logf_c2 sin logf_c, porque una curva cuadrática debe conservar también su término lineal. Para cada modelo se calcula R², R² ajustado, Cp de Mallows, AIC, BIC y PRESS.
Fórmulas: 8. Selección del modelo
Para un modelo con \(p\) parámetros (incluido \(\beta_0\)):
\[ R^2_p=1-\frac{SSE_p}{SST}, \qquad R^2_{aj,p}=1-\frac{n-1}{n-p}\cdot\frac{SSE_p}{SST} \]
\[ C_p=\frac{SSE_p}{\hat{\sigma}^2_{completo}}-(n-2p), \qquad \text{buen modelo si } C_p\approx p \]
\[ AIC_p=n\ln\!\left(\frac{SSE_p}{n}\right)+2p, \qquad BIC_p=n\ln\!\left(\frac{SSE_p}{n}\right)+p\ln n \]
\[ PRESS_p=\sum_{i=1}^{n}\left(\frac{e_i}{1-h_{ii}}\right)^2 \]
Criterio de elección: mayor \(R^2_{aj,p}\); \(C_p\) cercano a \(p\); menor \(AIC_p\), \(BIC_p\) y \(PRESS_p\).
Jerarquía: si el modelo incluye \((\log f_c)^2\), debe incluir también \(\log f_c\).
## 8. Seleccion del modelo: todas las regresiones posibles
# Estos tres terminos vienen de f y delta, que son positivos en NASA.xlsx.
# Centrar log(f) = restarle su media; el cuadrado permite curvatura.
logf_c <- log(f) - mean(log(f))
logf_c2 <- logf_c^2
logdelta <- log(delta)
datos2 <- data.frame(SSPL, logf_c, logf_c2, alpha, c, U_infinity, logdelta)
vars <- names(datos2)[-1]
full <- lm(SSPL ~ ., data = datos2)
SST <- sum((SSPL - mean(SSPL))^2) # suma de cuadrados total
MSE_k <- summary(full)$sigma^2 # CME del modelo completo
press <- function(m) sum((residuals(m) / (1 - hatvalues(m)))^2)
tabla <- NULL
for (r in 1:length(vars)) {
for (cmb in combn(vars, r, simplify = FALSE)) {
m <- lm(reformulate(cmb, "SSPL"), data = datos2)
p <- length(coef(m)) # parametros, incluido el intercepto
SSE <- sum(residuals(m)^2)
tabla <- rbind(tabla, data.frame(
variables = paste(cmb, collapse = " + "), p = p,
R2_p = 1 - SSE / SST, # R2
R2aj_p = 1 - ((n - 1) / (n - p)) * SSE / SST, # R2 ajustado
Cp = SSE / MSE_k - (n - 2 * p), # Cp de Mallows
AIC_p = n * log(SSE / n) + 2 * p,
BIC_p = n * log(SSE / n) + p * log(n),
PRESS = press(m)))
}
}
nrow(tabla) # 63 modelos[1] 63# Jerarquia: si esta logf_c2 debe estar logf_c.
tabla <- tabla[vapply(strsplit(tabla$variables, " + ", fixed = TRUE),
function(x) !("logf_c2" %in% x) || ("logf_c" %in% x),
logical(1)), ]
tabla <- tabla[order(tabla$BIC_p), ]
head(tabla, 10) variables p R2_p R2aj_p
60 logf_c + logf_c2 + c + U_infinity + logdelta 6 0.5368891 0.5353423
63 logf_c + logf_c2 + alpha + c + U_infinity + logdelta 7 0.5386730 0.5368227
57 logf_c + logf_c2 + alpha + c + U_infinity 6 0.5143835 0.5127616
59 logf_c + logf_c2 + alpha + U_infinity + logdelta 6 0.5018928 0.5002292
46 logf_c + logf_c2 + c + logdelta 5 0.4942210 0.4928705
58 logf_c + logf_c2 + alpha + c + logdelta 6 0.4942471 0.4925579
44 logf_c + logf_c2 + alpha + logdelta 5 0.4680206 0.4666001
51 logf_c + c + U_infinity + logdelta 5 0.4672206 0.4657980
61 logf_c + alpha + c + U_infinity + logdelta 6 0.4694556 0.4676835
47 logf_c + logf_c2 + U_infinity + logdelta 5 0.4597238 0.4582811
Cp AIC_p BIC_p PRESS
60 10.78467 4659.575 4691.466 33420.42
63 7.00000 4655.775 4692.981 33360.73
57 83.76628 4730.897 4762.788 35033.49
59 124.27128 4769.067 4800.958 35974.25
46 147.14965 4790.040 4816.616 36448.49
58 149.06501 4791.962 4823.853 36517.40
44 232.11276 4865.948 4892.525 38364.67
51 234.70696 4868.207 4894.783 38374.35
61 229.45947 4863.889 4895.780 38296.73
47 259.01801 4889.209 4915.785 38927.54tabla[which.max(tabla$R2aj_p), ] # mayor R2 ajustado variables p R2_p R2aj_p Cp
63 logf_c + logf_c2 + alpha + c + U_infinity + logdelta 7 0.538673 0.5368227 7
AIC_p BIC_p PRESS
63 4655.775 4692.981 33360.73tabla[which.min(abs(tabla$Cp - tabla$p)), ] # Cp mas cercano a p variables p R2_p R2aj_p Cp
63 logf_c + logf_c2 + alpha + c + U_infinity + logdelta 7 0.538673 0.5368227 7
AIC_p BIC_p PRESS
63 4655.775 4692.981 33360.73tabla[which.min(tabla$AIC_p), ] # menor AIC variables p R2_p R2aj_p Cp
63 logf_c + logf_c2 + alpha + c + U_infinity + logdelta 7 0.538673 0.5368227 7
AIC_p BIC_p PRESS
63 4655.775 4692.981 33360.73tabla[which.min(tabla$BIC_p), ] # menor BIC variables p R2_p R2aj_p Cp
60 logf_c + logf_c2 + c + U_infinity + logdelta 6 0.5368891 0.5353423 10.78467
AIC_p BIC_p PRESS
60 4659.575 4691.466 33420.42tabla[which.min(tabla$PRESS), ] # menor PRESS variables p R2_p R2aj_p Cp
63 logf_c + logf_c2 + alpha + c + U_infinity + logdelta 7 0.538673 0.5368227 7
AIC_p BIC_p PRESS
63 4655.775 4692.981 33360.73Interpretación de la tabla: R² y R² ajustado se prefieren altos; AIC, BIC y PRESS, bajos. En Cp interesa un valor cercano al número de parámetros p. PRESS suma los errores de predicción calculados al dejar fuera cada observación en el ajuste. AIC y BIC aquí omiten una constante común, por lo que se usan para comparar estos modelos en la misma base, no para igualar exactamente el número que devuelve AIC(m) o BIC(m).
Tras aplicar la jerarquía quedan 47 modelos válidos. El primero de la tabla ordenada tiene el BIC menor y usa: logf_c + logf_c2 + c + U_infinity + logdelta. Que un criterio prefiera un modelo no obliga a los demás a escoger el mismo.
9. Metodos paso a paso
step() hace una búsqueda más rápida que probar todas las combinaciones: forward empieza sin regresoras y las agrega; backward empieza con todas y las retira. Se muestran sus fórmulas finales como contraste con la búsqueda exhaustiva.
Fórmulas: 9. Métodos paso a paso
Forward: parte de \(\text{SSPL}=\beta_0+\varepsilon\) y agrega, en cada paso, la variable que más reduce el AIC. Backward: parte del modelo completo y quita, en cada paso, la variable cuya salida más reduce el AIC. Se detiene cuando ningún cambio lo reduce más.
## 9. Metodos paso a paso
# Usan las mismas seis columnas candidatas. El modelo final se elige de la
# tabla anterior, donde si se respeta la jerarquia del termino cuadratico.
nulo <- lm(SSPL ~ 1, data = datos2)
modelo_forward <- step(nulo, scope = formula(full),
direction = "forward", trace = 0)
modelo_backward <- step(full, direction = "backward", trace = 0)
formula(modelo_forward)SSPL ~ logf_c + logdelta + c + logf_c2 + U_infinity + alphaformula(modelo_backward)SSPL ~ logf_c + logf_c2 + alpha + c + U_infinity + logdeltaEstos procedimientos usan AIC y pueden llegar a modelos distintos. Además, step() no impone por sí solo la jerarquía de logf_c2; por eso la decisión final se toma de la tabla filtrada del punto 8.
10. Modelo elegido y conclusion
which.min(tabla$BIC_p) identifica el subconjunto válido con el BIC más bajo. Con esa fórmula se ajusta modelo_final. Sus supuestos se revisan nuevamente en el siguiente apartado.
Fórmulas: 10. Modelo elegido (menor BIC)
\[ \widehat{\text{SSPL}}=108.347-3.728\,\log f_c-1.445\,(\log f_c)^2-20.547\,c+0.0923\,U_\infty-3.153\,\log\delta \]
Para este modelo (\(k=5\)) se repiten las hipótesis y estadísticos de los apartados 2, 3, 4, 5, 6 y 7: el \(F_0\) con \(k=5\), el valor crítico \(\chi^2_{\alpha,5}\) de Breusch-Pagan, y Durbin-Watson y VIF con las fórmulas anteriores.
## 10. Modelo elegido y conclusion
# Menor BIC entre los subconjuntos validos de la seccion 8.
mejor <- tabla$variables[which.min(tabla$BIC_p)]
modelo_final <- lm(as.formula(paste("SSPL ~", mejor)), data = datos2)
summary(modelo_final)
Call:
lm(formula = as.formula(paste("SSPL ~", mejor)), data = datos2)
Residuals:
Min 1Q Median 3Q Max
-19.5973 -2.6729 -0.0638 3.0977 14.2376
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 108.347173 0.728268 148.77 <2e-16 ***
logf_c -3.727571 0.121081 -30.79 <2e-16 ***
logf_c2 -1.445308 0.096311 -15.01 <2e-16 ***
c -20.547492 1.301006 -15.79 <2e-16 ***
U_infinity 0.092342 0.007863 11.74 <2e-16 ***
logdelta -3.153082 0.107910 -29.22 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.703 on 1497 degrees of freedom
Multiple R-squared: 0.5369, Adjusted R-squared: 0.5353
F-statistic: 347.1 on 5 and 1497 DF, p-value: < 2.2e-16formula(modelo_final)SSPL ~ logf_c + logf_c2 + c + U_infinity + logdeltacat("Modelo final:", deparse(formula(modelo_final)), "\n")Modelo final: SSPL ~ logf_c + logf_c2 + c + U_infinity + logdelta cat("R2 ajustado:", summary(modelo_final)$adj.r.squared, "\n")R2 ajustado: 0.5353423 cat("Conclusiones: interpretar las asociaciones con cautela si fallan",
"los supuestos; el orden experimental no demuestra dependencia temporal.\n")Conclusiones: interpretar las asociaciones con cautela si fallan los supuestos; el orden experimental no demuestra dependencia temporal.El modelo seleccionado es SSPL ~ logf_c + logf_c2 + c + U_infinity + logdelta y su R² ajustado es 0.5353. Se escogió por BIC, no porque se haya demostrado que sus errores cumplen todos los supuestos.
Revisión de supuestos del modelo final
Se repite el mismo orden de revisión utilizado para el modelo original. Así puede evaluarse si la especificación seleccionada mejora los diagnósticos, sin asumir que una transformación los corrige automáticamente.
Normalidad
El Q-Q y el histograma muestran la forma de los residuales. Lilliefors, Anderson-Darling y Shapiro-Wilk se aplican ahora explícitamente al modelo final. La tabla compara el estadístico W y el valor p de Shapiro antes y después.
## Revision de supuestos del modelo final
## Supuesto de normalidad
# Verificacion grafica
res_final <- modelo_final$residuals
qqnorm(res_final)
qqline(res_final)
hist(res_final, breaks = 30, col = "steelblue",
main = "Histograma de residuales: modelo final")
# Verificacion formal
print(lillie.test(res_final))
Lilliefors (Kolmogorov-Smirnov) normality test
data: res_final
D = 0.040116, p-value = 6.63e-06print(ad.test(res_final))
Anderson-Darling normality test
data: res_final
A = 4.1511, p-value = 2.535e-10shapiro_original <- shapiro.test(residuales)
shapiro_final <- shapiro.test(res_final)
print(shapiro_final)
Shapiro-Wilk normality test
data: res_final
W = 0.98754, p-value = 4.595e-10tabla_shapiro <- data.frame(
Modelo = c("Original", "Final"),
W = c(unname(shapiro_original$statistic), unname(shapiro_final$statistic)),
Valor_p = c(shapiro_original$p.value, shapiro_final$p.value)
)
print(tabla_shapiro) Modelo W Valor_p
1 Original 0.9956327 2.465032e-04
2 Final 0.9875405 4.595129e-10
El estadístico de prueba de Lilliefors y su estadístico teórico para el modelo final se calculan igual que en el modelo original:
ord_f <- sort(res_final)
z_f <- (ord_f - mean(res_final)) / sd(res_final)
Fz_f <- pnorm(z_f)
i <- 1:n
D_final <- max(pmax(i / n - Fz_f, Fz_f - (i - 1) / n))
Dcrit_f <- 0.886 / sqrt(n)
data.frame(D_manual = D_final, D_teorico = Dcrit_f) D_manual D_teorico
1 0.04011566 0.02285358El modelo original obtuvo W = 0.9956 (p = 0.000247) y el final W = 0.9875 (p = 4.6e-10). Con α = 0.05, se rechaza la normalidad en el modelo final. Un W más próximo a 1 no basta por sí solo: también se revisan el valor p y el Q-Q.
Varianza constante
En la gráfica se busca una dispersión similar a lo largo de todos los valores ajustados. Se muestran Breusch-Pagan studentizado, la versión original y el valor crítico chi-cuadrado con los grados de libertad del modelo final, que pueden diferir de los cinco del modelo original.
## Supuesto de varianza constante
# Verificacion grafica
ajustados_final <- modelo_final$fitted.values
plot(x = ajustados_final, y = res_final)
abline(h = 0, col = "red")
# Verificacion formal
bp_final <- bptest(modelo_final)
cw_final <- bptest(modelo_final, studentize = FALSE)
print(bp_final)
studentized Breusch-Pagan test
data: modelo_final
BP = 331.4, df = 5, p-value < 2.2e-16print(cw_final)
Breusch-Pagan test
data: modelo_final
BP = 496.1, df = 5, p-value < 2.2e-16qchisq(0.05, length(coef(modelo_final)) - 1, lower.tail = FALSE)[1] 11.0705# Si BP supera el valor critico, se rechaza H0 de varianza constante.En la prueba studentizada, BP = 331.402 y p = <2e-16. Un p inferior a 0.05 aporta evidencia contra la varianza constante; una gráfica en forma de embudo también debe discutirse.
Independencia
El gráfico muestra si aparecen rachas por orden de observación. Durbin-Watson se presenta junto a su valor p; las cotas 1.718 y 1.820, citadas para unos 1500 registros y cinco regresoras, no se aplican automáticamente si el modelo final tiene otro número de términos. Además, el orden de ensayo no equivale a una serie temporal.
## Supuesto de independencia
# Verificacion grafica
orden <- 1:n
plot(x = orden, y = res_final)
abline(h = 0, col = "red")
# Verificacion formal
dw_final <- durbinWatsonTest(modelo_final, alternative = "two.sided")
print(dw_final) lag Autocorrelation D-W Statistic p-value
1 0.7879007 0.422193 0
Alternative hypothesis: rho != 0# Las cotas dl = 1.718 y du = 1.820 citadas corresponden aproximadamente
# a n = 1500, k = 5; si el modelo final cambia k, no se aplican sin otra tabla.Como esas cotas no se aplican sin otra tabla, se usa también la aproximación normal de \(d\):
d_final <- sum(diff(res_final)^2) / sum(res_final^2)
r1_final <- sum(res_final[-1] * res_final[-n]) / sum(res_final^2)
data.frame(d = d_final, r1 = r1_final,
limite_inferior = 2 - qnorm(0.975) * 2 / sqrt(n),
limite_superior = 2 + qnorm(0.975) * 2 / sqrt(n)) d r1 limite_inferior limite_superior
1 0.422193 0.7879007 1.898889 2.101111Multicolinealidad
Se revisa el VIF de los términos conservados. Un valor alto indica que la relación entre regresoras puede volver inestables los coeficientes; con una sola regresora, VIF no aporta una comparación entre variables.
## Multicolinealidad
if (length(coef(modelo_final)) > 2) vif(modelo_final) logf_c logf_c2 c U_infinity logdelta
1.171879 1.004477 1.005930 1.018346 1.161666 # VIF > 10: posible multicolinealidad elevada; con una regresora no aplica.Conclusión
A partir del análisis de la base de datos que usamos que fue NASA Airfoil Self-Noise, compuesta por 1.503 observaciones, se estudió la relación entre el nivel de presión sonora (SSPL) y cinco características aerodinámicas mediante un modelo de regresión lineal múltiple.
El modelo inicial fue estadísticamente significativo y explicó el 51,57 % de la variabilidad del nivel de presión sonora. Posteriormente, mediante la comparación de diferentes modelos con los criterios R² ajustado, Cp de Mallows, AIC, BIC y PRESS, se seleccionó el modelo con menor BIC, que alcanzó un R² de 0,5369 y un R² ajustado de 0,5353, mejorando ligeramente la capacidad explicativa del modelo inicial.
Sin embargo, la evaluación del modelo final evidenció problemas de normalidad y heterocedasticidad. También se identificó una posible autocorrelación en el orden de los registros, aunque debe interpretarse con precaución porque los datos no representan necesariamente una secuencia temporal. Por otra parte, los valores VIF no mostraron problemas importantes de multicolinealidad.
En conclusión, la regresión múltiple permitió identificar asociaciones entre las características aerodinámicas y el nivel de presión sonora. Aunque el modelo seleccionado presenta un ajuste estadísticamente significativo, se recomienda considerar errores estándar robustos, revisar las observaciones influyentes y estudiar la estructura de los datos. Por tanto, los resultados deben interpretarse con precaución debido a que no se cumplen todos los supuestos del modelo.
Referencias
Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723. https://doi.org/10.1109/TAC.1974.1100705
Allen, D. M. (1974). The relationship between variable selection and data augmentation and a method for prediction. Technometrics, 16(1), 125–127. https://doi.org/10.1080/00401706.1974.10489157
Anderson, T. W., & Darling, D. A. (1954). A test of goodness of fit. Journal of the American Statistical Association, 49(268), 765–769. https://doi.org/10.1080/01621459.1954.10501232
Breusch, T. S., & Pagan, A. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5), 1287–1294. https://doi.org/10.2307/1911963
Brooks, T., Pope, D., & Marcolini, M. (1989). Airfoil self-noise [Conjunto de datos]. UCI Machine Learning Repository. https://doi.org/10.24432/C5VW2C
Cook, R. D. (1977). Detection of influential observation in linear regression. Technometrics, 19(1), 15–18. https://doi.org/10.1080/00401706.1977.10489493
Durbin, J., & Watson, G. S. (1950). Testing for serial correlation in least squares regression. I. Biometrika, 37(3–4), 409–428. https://doi.org/10.2307/2332391
Durbin, J., & Watson, G. S. (1951). Testing for serial correlation in least squares regression. II. Biometrika, 38(1–2), 159–177. https://doi.org/10.2307/2332325
Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3.ª ed.). Sage.
Gross, J., & Ligges, U. (2015). nortest: Tests for normality (Paquete de R, versión 1.0-4). https://CRAN.R-project.org/package=nortest
Koenker, R. (1981). A note on studentizing a test for heteroscedasticity. Journal of Econometrics, 17(1), 107–112. https://doi.org/10.1016/0304-4076(81)90062-2
Lilliefors, H. W. (1967). On the Kolmogorov-Smirnov test for normality with mean and variance unknown. Journal of the American Statistical Association, 62(318), 399–402. https://doi.org/10.1080/01621459.1967.10482916
Mallows, C. L. (1973). Some comments on Cp. Technometrics, 15(4), 661–675. https://doi.org/10.1080/00401706.1973.10489103
Montgomery, D. C., Peck, E. A., & Vining, G. G. (2021). Introduction to linear regression analysis (6.ª ed.). Wiley.
R Core Team. (s. f.). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/
Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464. https://doi.org/10.1214/aos/1176344136
Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3–4), 591–611. https://doi.org/10.1093/biomet/52.3-4.591
Wickham, H., & Bryan, J. (s. f.). readxl: Read Excel files [Paquete de R]. https://CRAN.R-project.org/package=readxl
Zeileis, A., & Hothorn, T. (2002). Diagnostic checking in regression relationships. R News, 2(3), 7–10.
Zeileis, A., & Hothorn, T. (2002). Diagnostic checking in regression relationships. R News, 2(3), 7–10.