# 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 <- 5Trabajo final — Regresión lineal múltiple
NASA Airfoil Self-Noise
Introducción
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 permite ir directamente a cada apartado. Todos los bloques de código se muestran en el informe.
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.
Exploracion
Se revisan resumen, valores faltantes, correlaciones y relaciones gráficas entre variables.
## 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)
1. Modelo y estimacion de parametros
Se estima por mínimos cuadrados el modelo lineal completo con las cinco regresoras originales.
## 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-162. Inferencia sobre el modelo completo
Se contrasta si todas las pendientes son cero frente a que al menos una es distinta de cero.
## 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.2200763. Inferencia para cada parametro
Se examinan los estadísticos t y los valores p de cada coeficiente.
## 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+024. Normalidad de los residuales
El gráfico Q-Q y las pruebas formales examinan la distribución de los errores.
## 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.00024655. Varianza constante
El gráfico de residuos y la prueba de Breusch-Pagan evalúan la homocedasticidad.
## 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.07056. Independencia
Se observa el orden de las mediciones y se aplica Durbin-Watson con cautela: el orden experimental no es una serie temporal.
## 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.7. Multicolinealidad
El VIF mide cuánto se relaciona cada regresora con las demás.
## 7. Multicolinealidad
vif(modelo) f alpha c U_infinity delta
1.144444 3.441658 1.510754 1.041698 2.532127 7.1. Observaciones influyentes
La distancia de Cook y los residuos estandarizados señalan observaciones que conviene revisar.
## 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] 138. Seleccion del modelo: todas las regresiones posibles
Se forman logf_c = log(f) - media(log(f)), logf_c2 = logf_c² y logdelta = log(delta). Son transformaciones de variables medidas, no datos nuevos. Se comparan 63 subconjuntos; luego se excluyen los que incluyen el cuadrado sin el término lineal. R², R² ajustado, Cp, AIC, BIC y PRESS permiten comparar los modelos.
## 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.739. Metodos paso a paso
Se muestran selección hacia adelante y eliminación hacia atrás. El modelo final se decide con la tabla jerárquica del punto 8.
## 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 + logdelta10. Modelo elegido y conclusion
Se selecciona el menor BIC entre los subconjuntos válidos y se repiten los diagnósticos principales. Un BIC menor no garantiza que los supuestos se cumplan.
## 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 + logdelta# El modelo final debe volver a revisarse: seleccionar variables no garantiza
# normalidad, homocedasticidad ni independencia de los errores.
res_final <- residuals(modelo_final)
qqnorm(res_final, main = "Q-Q del modelo final")
qqline(res_final, col = "red")
plot(fitted(modelo_final), res_final, pch = 20,
main = "Residuales vs ajustados: modelo final",
xlab = "Valores ajustados", ylab = "Residuales")
abline(h = 0, col = "red")
print(bptest(modelo_final))
studentized Breusch-Pagan test
data: modelo_final
BP = 331.4, df = 5, p-value < 2.2e-16print(durbinWatsonTest(modelo_final, alternative = "two.sided")) lag Autocorrelation D-W Statistic p-value
1 0.7879007 0.422193 0
Alternative hypothesis: rho != 0if (length(coef(modelo_final)) > 2) print(vif(modelo_final)) logf_c logf_c2 c U_infinity logdelta
1.171879 1.004477 1.005930 1.018346 1.161666 cat("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.Cierre
Interpreta los coeficientes y la selección únicamente después de revisar las salidas que genera el código. Si los diagnósticos señalan problemas, presenta las limitaciones antes de sacar conclusiones sobre la relación entre las variables. El informe es autocontenido una vez tejido en HTML.