Informe 2 - Regresión lineal multiple
Jane Juliana Gutiérrez Londoño (Grupo 2), Nicolas Valencia Trujillo (Grupo 2), Pablo Restrepo Osorio (Grupo 2)
04 de septiembre de 2021
require(car)## Loading required package: car## Warning: package 'car' was built under R version 4.0.5## Loading required package: carDatalibrary(perturb)
require(leaps)## Loading required package: leaps## Warning: package 'leaps' was built under R version 4.0.5library(corrplot)## corrplot 0.90 loadedlibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.0.4##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondatos <- read.csv("datos2.csv", sep = " ")\(1.\) Estime un modelo de regresión lineal múltiple que explique el Riesgo de Infección en términos de todas las variables predictoras. Analice la significancia de la regresión y de los parámetros individuales. Interprete los parámetros estimados. Calcule e interprete el coeficiente de determinación múltiple R2.Comente los resultados.
miAnova=function(modeloreg){
SSq=unlist(anova(modeloreg)["Sum Sq"])
k=length(SSq)-1
SSR=sum(SSq[1:k])
SSE=SSq[(k+1)]
MSR=SSR/k
df.error=unlist(anova(modeloreg)["Df"])[k+1]
MSE=SSE/df.error
F0=MSR/MSE
VP=pf(F0,k,df.error,lower.tail=F)
result=data.frame(SumSq=c(SSR,SSE),Df=c(k,df.error),MeanSq=c(MSR,MSE),F0=c(round(F0,digits=3),' '),
P.value=c(format(VP,scientific =TRUE,digits=3),' '),row.names=c("Modelo","Error"))
cat("Tabla ANOVA Modelo de Regresión","\n")
result}attach(datos)
modelo<-lm(Y~X1+X2+X3+X4+X5)coeficientes de nuestro modelo de regresión
summary(modelo)##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.03944 -0.80043 -0.00266 0.60450 2.23292
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5986009 1.5159559 -0.395 0.694365
## X1 0.2106683 0.0785765 2.681 0.009501 **
## X2 0.0197512 0.0277108 0.713 0.478803
## X3 0.0470925 0.0132888 3.544 0.000779 ***
## X4 0.0105604 0.0073166 1.443 0.154213
## X5 0.0008996 0.0007379 1.219 0.227679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.035 on 59 degrees of freedom
## Multiple R-squared: 0.4468, Adjusted R-squared: 0.3999
## F-statistic: 9.53 on 5 and 59 DF, p-value: 1.058e-06Ahora entonces nuestro modelo estimado sería:
\[{Y}_{}=-0.599+0.211X_{1}+0.02X_{2}+0.047X_{3}+0.011X_{4}+0.001X_{5}\]
Siendo \(Y\) el Riesgo de infección, \(X_{1}\) la Duración de la estadía, \(X_{2}\) la Rutina de cultivos, \(X_{3}\) el Número de camas, \(X_{4}\) el Censo promedio diario, \(X_{5}\) el Número de enfermeras.
Para analizar la significancia de la regresión \(\alpha =0.05\) y el valor p ahora, planteemos nuestra prueba de significancia y veamos la tabla.
Hipótesis: \(H_{0}:\beta _{1}=\beta _{2}=\beta _{3}=\beta _{4}=\beta _{5}= 0\) vs \(H_{1}:\)Al menos un\(\beta _{i}\neq 0\) para \(\beta _{i}\) con i=0,1,2,3,4,5.
Estadístico de prueba: \(F_{0}=\frac{MSR}{MSE} \sim ^{Bajo H_{0}} f_{0.05,5,65-6}\)
Se rechaza si el \(valor\) \(p=P(f_{0.05,5,65-6}>F_{0})\) es pequeño.
miAnova(modelo)## Tabla ANOVA Modelo de Regresión## SumSq Df MeanSq F0 P.value
## Modelo 51.04527 5 10.209053 9.53 1.06e-06
## Error 63.20335 59 1.071243Para nuestra regresión podemos ver que el valor p es muy pequeño, por lo tanto nuestra regresión es significativa.
Para analizar la significancia de los parametros individuales utilizaremos un \[\alpha =0.05\] y el valor p ahora, planteemos nuestra prueba de significancia y utilizemos los valores calculados en la primera tabla.
Para cada \(\beta _{i}\) con i=0,1,2,3,4,5:
Hipótesis: \(H_{0}:\beta _{i}= 0\) vs \(H_{1}:\beta _{i}\neq 0\)
Estadístico de prueba: \(t_{0}=\frac{\hat{\beta _{i}}}{Se(\beta _{i})} \sim ^{Bajo H_{0}} t_{65-6}\)
Se rechaza si el \(valor p=P(t_{\frac{0.05 }{2},65-6}>|t_{0}|)\) es pequeño.
El valor p de \(\beta _{0}\) no pequeno esta no es sinificatica y no se puede interpretar.
El valor p de \(\beta _{1}\) es pequeño entonces es significativa podemos interpretarla como que pir un aumento unitario en X1 se espera que Y aumente en promedio 0.211 unidades, siempre que las demás variables predictoras permanezcan constantes.
El valor p de \(\beta _{2}\) no pequeno esta no es sinificatica y no se puede interpretar.
El valor p de \(\beta _{3}\) es pequeño entonces es significativa podemos interpretarla como que pir un aumento unitario en X1 se espera que Y aumente en promedio 0.047 unidades, siempre que las demás variables predictoras permanezcan constantes.
El valor p de \(\beta _{4}\) no pequeno esta no es sinificatica y no se puede interpretar.
El valor p de \(\beta _{5}\) no pequeno esta no es sinificatica y no se puede interpretar.
Y por ultimo el valor del coeficiente de determinación multiple \(R^{2}\) es 0.4468, lo que quiere decir que un 44,68% de la variabilidad de la variable respuesta es explicada por la regresión.
\(2.\) Use la tabla de todas las regresiones posibles, para probar la significancia simultánea del subconjunto de tres variables con los valores p mayores del punto anterior. Según el resultado de la prueba es posible descartar del modelo las variables del subconjunto?.
attach(datos)## The following objects are masked from datos (pos = 3):
##
## X1, X2, X3, X4, X5, Ymodeloreg = lm(Y~X1+X2+X3+X4+X5)summary(modeloreg)##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.03944 -0.80043 -0.00266 0.60450 2.23292
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5986009 1.5159559 -0.395 0.694365
## X1 0.2106683 0.0785765 2.681 0.009501 **
## X2 0.0197512 0.0277108 0.713 0.478803
## X3 0.0470925 0.0132888 3.544 0.000779 ***
## X4 0.0105604 0.0073166 1.443 0.154213
## X5 0.0008996 0.0007379 1.219 0.227679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.035 on 59 degrees of freedom
## Multiple R-squared: 0.4468, Adjusted R-squared: 0.3999
## F-statistic: 9.53 on 5 and 59 DF, p-value: 1.058e-06myAllRegTable <- function(lm.model, response = model.response(model.frame(lm.model)), MSE = F){
regTable <- summary(regsubsets(model.matrix(lm.model)[, -1], response,
nbest = 2^(lm.model$rank - 1) - 1, really.big = T))
pvCount <- as.vector(apply(regTable$which[, -1], 1, sum))
pvIDs <- apply(regTable$which[, -1], 1, function(x) as.character(paste(colnames(model.matrix(lm.model)[, -1])[x],
collapse = " ")))
result <- if(MSE){
data.frame(k = pvCount, R_sq = round(regTable$rsq, 3), adj_R_sq = round(regTable$adjr2, 3),
MSE = round(regTable$rss/(nrow(model.matrix(lm.model)[,-1]) - (pvCount + 1)), 3),
Cp = round(regTable$cp, 3), Variables_in_model = pvIDs)
} else {
data.frame(k = pvCount, R_sq = round(regTable$rsq, 3), adj_R_sq = round(regTable$adjr2, 3),
SSE = round(regTable$rss, 3),
Cp = round(regTable$cp, 3), Variables_in_model = pvIDs)
}
format(result, digits = 6)
}myAllRegTable(modeloreg)## k R_sq adj_R_sq SSE Cp Variables_in_model
## 1 1 0.249 0.237 85.813 19.106 X3
## 2 1 0.244 0.232 86.376 19.632 X1
## 3 1 0.182 0.169 93.468 26.252 X4
## 4 1 0.068 0.053 106.508 38.425 X5
## 5 1 0.001 -0.014 114.090 45.503 X2
## 6 2 0.411 0.392 67.257 3.784 X1 X3
## 7 2 0.319 0.297 77.804 13.629 X3 X4
## 8 2 0.311 0.289 78.739 14.503 X1 X4
## 9 2 0.293 0.270 80.764 16.392 X3 X5
## 10 2 0.276 0.253 82.667 18.169 X2 X3
## 11 2 0.258 0.235 84.721 20.087 X1 X5
## 12 2 0.248 0.224 85.877 21.166 X1 X2
## 13 2 0.234 0.209 87.519 22.698 X4 X5
## 14 2 0.182 0.156 93.416 28.203 X2 X4
## 15 2 0.071 0.041 106.102 40.046 X2 X5
## 16 3 0.430 0.402 65.108 3.778 X1 X3 X4
## 17 3 0.422 0.393 66.088 4.693 X1 X3 X5
## 18 3 0.415 0.386 66.850 5.404 X1 X2 X3
## 19 3 0.359 0.327 73.279 11.405 X3 X4 X5
## 20 3 0.336 0.303 75.868 13.822 X2 X3 X4
## 21 3 0.328 0.295 76.791 14.684 X1 X4 X5
## 22 3 0.325 0.292 77.090 14.963 X2 X3 X5
## 23 3 0.314 0.280 78.400 16.186 X1 X2 X4
## 24 3 0.261 0.224 84.459 21.842 X1 X2 X5
## 25 3 0.236 0.198 87.321 24.514 X2 X4 X5
## 26 4 0.442 0.405 63.748 4.508 X1 X3 X4 X5
## 27 4 0.433 0.395 64.795 5.486 X1 X2 X3 X4
## 28 4 0.427 0.389 65.435 6.083 X1 X2 X3 X5
## 29 4 0.379 0.338 70.904 11.188 X2 X3 X4 X5
## 30 4 0.329 0.284 76.656 16.558 X1 X2 X4 X5
## 31 5 0.447 0.400 63.203 6.000 X1 X2 X3 X4 X5Prueba de la significancia del subconjunto formado por X2,X4 y X5
La significancia simultánea de estas variables la vamos a probar por medio de una hipótesis:
\(H_{0}:\) \(\beta(X_{2}+X_{4}+X_{5}) , \beta = 0\)
vs \(H_{1}:\) \(\beta_{2}X_{2}+\beta_{4}X_{4}+\beta_{5}X_{5}, \beta_{2} \mid \beta_{4} \mid \beta_{5} \neq 0\)
\(M.C. = Y_{i}=\beta_{0}+\beta_{1}X_{i1}+\beta_{2}X_{i2}+\beta_{3}X_{i3}+\beta_{4}X_{i4}+\beta_{5}X_{i5}+\varepsilon_{i}\)
\(M.R. = Y_{i}=\beta_{0}+\beta_{1}X_{i1}+\beta_{3}X_{i3}+\varepsilon_{i}\)
Estadístico de prueba: \(\displaystyle F_{0} = \frac{\frac{SSE(M.R.)-SSE(M.C.)}{r}}{MSE(M.C.)} \sim f_{r,n-p}\)
\(SSE(M.R.) = 67.257\)
\(SSE(M.C.) = 63.203\)
\(MSE(M.C.) = \frac{63.203}{65-6}\)
\(r = 3\)
\(n = 65\)
\(p = 6\)
\(F_{0} = \frac{\frac{67.257-63.203}{3}}{\frac{63.203}{65-6}} = 1.261\)
Se rechaza \(H_{0}\) si \(F_{0}>f_{0.05,3,59}\)
\(f_{0.05,3,59} = 0.1167\)
\(1.261>0.1167\)
Por ende se rechaza \(H_{0}\) y no se puede descartar este subconjunto del modelo enteramente.
\(3.\) Plantee una pregunta donde su solución implique el uso exclusivo de una prueba de hipótesis lineal general de la forma H0 : Lβ = 0 (solo se puede usar este procedimiento y no SSextra), donde especifique claramente la matriz L, el modelo reducido y la expresión para el estadístico de prueba.
Hipótesis Lineal General
Queremos probar la siguiente hipótesis:
- \(H_{0}: \beta_{1} = 0 \;,\; \beta_{2} = \beta_{3}\;,\; \beta_{4} = \beta_{5}\)
vs
- \(H_{1}: \beta_{1} \neq 0 \;o\; \beta_{2} \neq \beta_{3}\;o\; \beta_{4} \neq \beta_{5}\)
Igualamos las ecuaciones a cero y tenemos:
\(H_{0}: \beta_{1} = 0 \;,\; \beta_{2}-\beta_{3} = 0 \;,\; \beta_{4} - \beta_{5} = 0\)
\(H_{1}: \beta_{1} \neq 0 \;o\; \beta_{2} - \beta_{3} \neq 0\;o\; \beta_{4} - \beta_{5} \neq 0\)
Armamos la matriz L:
\(L = \begin{pmatrix} 0 &1 &0 &0 &0 &0 \\ 0&0 &1 &-1 &0 &0 \\ 0 & 0 & 0 &0 &1 &-1 \end{pmatrix}\)
\(M.C: = Y_{i}=\beta_{0}+\beta_{1}X_{i1}+\beta_{2}X_{i2}+\beta_{3}X_{i3}+\beta_{4}X_{i4}+\beta_{5}X_{i5}+\varepsilon_{i}\)
\(M.R: Y_{i}=\beta_{0}+\beta_{2}Z_{i2,3}+\beta_{4}Z_{i4,5} + \varepsilon_{i}\)
Estadístico de prueba: \(\displaystyle F_{0} = \frac{\frac{SSE(M.R.)-SSE(M.C.)}{m}}{MSE(M.C.)} \sim f_{m,n-p}\)
\(4.\) Realice una validación de los supuestos en los errores y examine si hay valores atípicos, de balanceo e influenciales. Qué puede decir acerca de la validez de éste modelo?. Argumente.
Validación de supuestos del modelo
Supuesto de Media Cero
ei = modelo$residuals
yi_mod = modelo$fitted.values
round(mean(ei),3)## [1] 0Gracias a esto validamos que la media es efectivamente 0, se valida el supuesto de Media Cero.
Supuesto de Varianza Constante
plot(yi_mod,ei,main="Residuales vs. Valores Ajustados", xlab = "Valores Ajustados", ylab = "Residuales")
abline(h=0,lty=2,col=2)
require(ggplot2)## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.0.4
data_val=data.frame(yi_mod,ei)
ggplot(data_val,aes(x=yi_mod,y=ei))+geom_point()+theme_bw()+geom_smooth()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Aquí podemos observar que no se estaría cumpliendo gráficamente los requerimientos para validar una varianza constante, el gráfico sólo nos dice que sospechemos la no linealidad del modelo, no se valida el supuesto de Varianza constante.
Supuesto de Normalidad
ei## 1 2 3 4 5 6
## -0.41425353 -1.51713870 0.24297542 1.72677333 0.92509934 0.35521408
## 7 8 9 10 11 12
## -1.60590283 -0.05386955 0.26767225 0.29864785 0.42082173 0.60450491
## 13 14 15 16 17 18
## 1.58629960 0.12311898 -0.64599426 0.76809511 -0.84252589 0.04781284
## 19 20 21 22 23 24
## -1.17309945 -1.02080307 0.93338045 -0.68622392 0.49346944 0.44652340
## 25 26 27 28 29 30
## -0.86352827 0.38108073 -0.49114604 -0.78238075 1.19336671 -1.20564532
## 31 32 33 34 35 36
## 0.59383601 -0.92949453 0.57338729 -0.00265947 -0.92053017 -0.37771269
## 37 38 39 40 41 42
## -2.03943547 -0.44898491 0.60874392 -0.80042500 1.61161226 2.23291925
## 43 44 45 46 47 48
## 1.58161195 -0.26441992 1.27181927 1.11795399 0.22967564 -0.24346947
## 49 50 51 52 53 54
## -0.20274140 -1.03970211 1.50474011 -1.59021223 0.02563588 -0.74234671
## 55 56 57 58 59 60
## -1.53179487 0.87718547 -1.30584847 -0.21030709 -0.36210096 2.09343412
## 61 62 63 64 65
## -0.19630494 -1.13707927 0.16500123 -0.83717797 1.18284669hist(ei,col="gray", main = "Histograma de Residuales", xlab ="Residuales", ylab = "Frecuencia")
qqnorm(ei)
qqline(ei,col="red")
shapiro.test(ei)##
## Shapiro-Wilk normality test
##
## data: ei
## W = 0.98462, p-value = 0.5979alpha = 0.05Gráfica
No hay evidencia fuerte en contra del supuesto de que los residuales distribuyen normal.
Shapiro Wilk
0.5979. No rechaza H0 que dice que se cumple el supuesto de normalidad.
Estadísticos de salida
t1<-predict(modelo,se.fit=T)
t2<-round(residuals(modelo),4)
t3<-round(cooks.distance(modelo),4)
t4<-round(hatvalues(modelo),4)
t5<-round(dffits(modelo),4)
restud<-round(rstudent(modelo),4)
est_salida <- data.frame(datos$Y,yhat=round(t1$fit,4),
se.yhat=round(t1$se.fit,6),residuals=t2
,res.estud=restud,Cooks.D=t3,hii.value=t4,Dffits=t5)
(est_salida)## datos.Y yhat se.yhat residuals res.estud Cooks.D hii.value Dffits
## 1 3.7 4.1143 0.281627 -0.4143 -0.4130 0.0023 0.0740 -0.1168
## 2 2.8 4.3171 0.223854 -1.5171 -1.5179 0.0184 0.0468 -0.3362
## 3 4.2 3.9570 0.193758 0.2430 0.2371 0.0003 0.0350 0.0452
## 4 6.2 4.4732 0.306929 1.7268 1.7787 0.0490 0.0879 0.5523
## 5 5.7 4.7749 0.352663 0.9251 0.9499 0.0198 0.1161 0.3443
## 6 4.5 4.1448 0.265863 0.3552 0.3525 0.0015 0.0660 0.0937
## 7 1.6 3.2059 0.293548 -1.6059 -1.6411 0.0382 0.0804 -0.4854
## 8 5.1 5.1539 0.287489 -0.0539 -0.0537 0.0000 0.0772 -0.0155
## 9 4.1 3.8323 0.218296 0.2677 0.2625 0.0005 0.0445 0.0566
## 10 4.4 4.1014 0.237148 0.2986 0.2941 0.0008 0.0525 0.0692
## 11 5.0 4.5792 0.162106 0.4208 0.4088 0.0007 0.0245 0.0648
## 12 4.3 3.6955 0.249986 0.6045 0.5986 0.0037 0.0583 0.1490
## 13 5.3 3.7137 0.201230 1.5863 1.5822 0.0160 0.0378 0.3136
## 14 4.8 4.6769 0.402158 0.1231 0.1280 0.0005 0.1510 0.0540
## 15 4.4 5.0460 0.185804 -0.6460 -0.6312 0.0022 0.0322 -0.1152
## 16 5.3 4.5319 0.318533 0.7681 0.7774 0.0106 0.0947 0.2514
## 17 2.9 3.7425 0.213885 -0.8425 -0.8298 0.0051 0.0427 -0.1753
## 18 4.3 4.2522 0.356599 0.0478 0.0488 0.0001 0.1187 0.0179
## 19 2.0 3.1731 0.259816 -1.1731 -1.1747 0.0154 0.0630 -0.3046
## 20 2.7 3.7208 0.314166 -1.0208 -1.0358 0.0181 0.0921 -0.3300
## 21 5.6 4.6666 0.343458 0.9334 0.9553 0.0188 0.1101 0.3360
## 22 4.1 4.7862 0.463640 -0.6862 -0.7387 0.0230 0.2007 -0.3701
## 23 6.6 6.1065 0.490120 0.4935 0.5380 0.0141 0.2242 0.2893
## 24 5.1 4.6535 0.189494 0.4465 0.4358 0.0011 0.0335 0.0812
## 25 4.5 5.3635 0.287627 -0.8635 -0.8667 0.0105 0.0772 -0.2507
## 26 4.3 3.9189 0.398923 0.3811 0.3962 0.0046 0.1486 0.1655
## 27 6.5 6.9911 0.692102 -0.4911 -0.6350 0.0549 0.4471 -0.5711
## 28 2.9 3.6824 0.421061 -0.7824 -0.8252 0.0226 0.1655 -0.3675
## 29 4.5 3.3066 0.280221 1.1934 1.2022 0.0189 0.0733 0.3381
## 30 4.9 6.1056 0.464484 -1.2056 -1.3114 0.0714 0.2014 -0.6586
## 31 5.6 5.0062 0.213415 0.5938 0.5831 0.0025 0.0425 0.1229
## 32 3.0 3.9295 0.380628 -0.9295 -0.9652 0.0243 0.1352 -0.3817
## 33 5.7 5.1266 0.381952 0.5734 0.5928 0.0093 0.1362 0.2354
## 34 5.0 5.0027 0.208808 -0.0027 -0.0026 0.0000 0.0407 -0.0005
## 35 2.9 3.8205 0.230729 -0.9205 -0.9110 0.0073 0.0497 -0.2083
## 36 4.5 4.8777 0.208361 -0.3777 -0.3698 0.0010 0.0405 -0.0760
## 37 2.5 4.5394 0.270130 -2.0394 -2.0993 0.0508 0.0681 -0.5676
## 38 3.4 3.8490 0.271107 -0.4490 -0.4464 0.0025 0.0686 -0.1212
## 39 5.8 5.1913 0.403037 0.6087 0.6353 0.0121 0.1516 0.2686
## 40 4.8 5.6004 0.307301 -0.8004 -0.8075 0.0106 0.0882 -0.2511
## 41 5.4 3.7884 0.409678 1.6116 1.7237 0.0890 0.1567 0.7429
## 42 6.3 4.0671 0.146814 2.2329 2.2535 0.0163 0.0201 0.3229
## 43 6.3 4.7184 0.291469 1.5816 1.6141 0.0364 0.0793 0.4737
## 44 3.4 3.6644 0.418980 -0.2644 -0.2772 0.0025 0.1639 -0.1227
## 45 7.8 6.5282 0.501587 1.2718 1.4167 0.1010 0.2349 0.7849
## 46 6.4 5.2820 0.234490 1.1180 1.1112 0.0111 0.0513 0.2585
## 47 4.6 4.3703 0.191635 0.2297 0.2240 0.0003 0.0343 0.0422
## 48 3.1 3.3435 0.241952 -0.2435 -0.2400 0.0006 0.0546 -0.0577
## 49 4.1 4.3027 0.156387 -0.2027 -0.1965 0.0002 0.0228 -0.0300
## 50 2.9 3.9397 0.285121 -1.0397 -1.0458 0.0149 0.0759 -0.2997
## 51 4.7 3.1953 0.276387 1.5047 1.5255 0.0291 0.0713 0.4227
## 52 5.4 6.9902 0.606346 -1.5902 -1.9397 0.3130 0.3432 -1.4021
## 53 4.8 4.7744 0.240466 0.0256 0.0252 0.0000 0.0540 0.0060
## 54 2.3 3.0423 0.246426 -0.7423 -0.7356 0.0055 0.0567 -0.1803
## 55 2.0 3.5318 0.209103 -1.5318 -1.5281 0.0162 0.0408 -0.3152
## 56 5.5 4.6228 0.252506 0.8772 0.8721 0.0081 0.0595 0.2194
## 57 1.4 2.7058 0.298558 -1.3058 -1.3261 0.0263 0.0832 -0.3995
## 58 4.7 4.9103 0.201498 -0.2103 -0.2055 0.0003 0.0379 -0.0408
## 59 3.9 4.2621 0.211224 -0.3621 -0.3547 0.0009 0.0416 -0.0739
## 60 5.5 3.4066 0.215746 2.0934 2.1291 0.0324 0.0435 0.4538
## 61 3.7 3.8963 0.213522 -0.1963 -0.1922 0.0003 0.0426 -0.0405
## 62 3.9 5.0371 0.224410 -1.1371 -1.1280 0.0104 0.0470 -0.2505
## 63 4.2 4.0350 0.378101 0.1650 0.1698 0.0008 0.1335 0.0667
## 64 2.9 3.7372 0.372875 -0.8372 -0.8652 0.0187 0.1298 -0.3342
## 65 5.6 4.4172 0.156936 1.1828 1.1596 0.0052 0.0230 0.1779Valores atípicos.
Bajo el criterio \(|r_𝒊|>3\)
library(dplyr)
observaciones_atipicas <- filter(est_salida, est_salida$res.estud > 3 | est_salida$res.estud < -3)
observaciones_atipicas## [1] datos.Y yhat se.yhat residuals res.estud Cooks.D hii.value
## [8] Dffits
## <0 rows> (or 0-length row.names)No hay ningun residual estudentizado que en valor absoluto sea mayor a 3, por tanto podemos concluir que no hay observaciones atipicas en nuestro modelo.
Puntos de balanceo.
Bajo el criterio \(ℎ_{ii} > \frac{2𝑝}{𝑛}\)
puntos_balanceo <- filter(est_salida, as.double(est_salida$hii.value) > (2*6 / 65))
puntos_balanceo## datos.Y yhat se.yhat residuals res.estud Cooks.D hii.value Dffits
## 22 4.1 4.7862 0.463640 -0.6862 -0.7387 0.0230 0.2007 -0.3701
## 23 6.6 6.1065 0.490120 0.4935 0.5380 0.0141 0.2242 0.2893
## 27 6.5 6.9911 0.692102 -0.4911 -0.6350 0.0549 0.4471 -0.5711
## 30 4.9 6.1056 0.464484 -1.2056 -1.3114 0.0714 0.2014 -0.6586
## 45 7.8 6.5282 0.501587 1.2718 1.4167 0.1010 0.2349 0.7849
## 52 5.4 6.9902 0.606346 -1.5902 -1.9397 0.3130 0.3432 -1.4021Se identifico que las observaciones 22, 23, 27, 30, 45 y 52 son puntos de balanceo.
Puntos influyentes
Bajo el criterio \(|𝐷𝐹𝐹𝐼𝑇𝑆 _𝑖| > 2\sqrt\frac{𝑝}{𝑛}\)
puntos_influyentes <- filter(est_salida, est_salida$Dffits > 2*sqrt(6/65) | est_salida$Dffits < -sqrt(2*6/65))
puntos_influyentes## datos.Y yhat se.yhat residuals res.estud Cooks.D hii.value Dffits
## 7 1.6 3.2059 0.293548 -1.6059 -1.6411 0.0382 0.0804 -0.4854
## 27 6.5 6.9911 0.692102 -0.4911 -0.6350 0.0549 0.4471 -0.5711
## 30 4.9 6.1056 0.464484 -1.2056 -1.3114 0.0714 0.2014 -0.6586
## 37 2.5 4.5394 0.270130 -2.0394 -2.0993 0.0508 0.0681 -0.5676
## 41 5.4 3.7884 0.409678 1.6116 1.7237 0.0890 0.1567 0.7429
## 45 7.8 6.5282 0.501587 1.2718 1.4167 0.1010 0.2349 0.7849
## 52 5.4 6.9902 0.606346 -1.5902 -1.9397 0.3130 0.3432 -1.4021Bajo este criterio encontramos que las observaciones 7, 27, 30, 37, 41, 45 y 52 son puntos influyentes.
Bajo el criterio \(cooks.di > 1\)
puntos_influyentes2 <- filter(est_salida, est_salida$Cooks.D > 1)
puntos_influyentes2## [1] datos.Y yhat se.yhat residuals res.estud Cooks.D hii.value
## [8] Dffits
## <0 rows> (or 0-length row.names)No se encontro ningún punto influyente con las distancias de cook
Validez del modelo
El modelo cumple con el supuesto de media cero y no hay fuerte evidencia en contra del supuesto de normalidad. El único supuesto que no se cumple es el de varianza constante y hay sospecha de no linealidad por tanto el modelo no es totalmente valido.
\(5.\) Verificar la presencia de multicolinealidad usando graficos y/o indicadores apropiados.
panel.cor=function(x,y,digits=2,prefix="",cex.cor){
usr=par("usr");on.exit(par(usr))
par(usr=c(0,1,0,1))
r=cor(x,y)
txt=format(c(r,0.123456789),digits=digits)[1]
txt=paste(prefix,txt,sep="")
if(missing(cex.cor))
cex = 0.4/strwidth(txt)
text(0.5, 0.5, txt, cex = cex)}miscoeficientes=function(modeloreg,datosreg){
coefi=coef(modeloreg)
datos=as.data.frame(scale(datosreg))
coef.std=c(0,coef(lm(update(formula(modeloreg),~.+0),datos)))
limites=confint(modeloreg,level=0.95)
vifs=c(0,vif(modeloreg))
resul=data.frame(Estimación=coefi,Limites=limites,Vif=vifs,Coef.Std=coef.std)
cat("Coeficientes estimados, sus I.C, Vifs y Coeficientes estimados estandarizados","\n")
resul}misDiagnostcolin=function(modeloreg,centrar=F){
if(centrar==F){
X=model.matrix(modeloreg)
val.prop=prcomp(X,center=FALSE,scale=TRUE)$sdev^2
Ind=colldiag(modeloreg)
resul=data.frame(Val.propio=val.prop,Ind.Cond=Ind$condindx,Pi=Ind$pi)
cat("Diagnósticos Multicolinealidad - Intercepto incluído","\n", "Índices de Condición y Proporciones deVarianza","\n")
}else{
X=model.matrix(modeloreg)[,-1]
val.prop=prcomp(X,center=TRUE,scale=TRUE)$sdev^2
Ind=colldiag(modeloreg,center=TRUE,scale=TRUE)
resul=data.frame(Val.propio=val.prop,Ind.Cond=Ind$condindx,Pi=Ind$pi)
cat("Diagnósticos Multicolinealidad -Intercepto ajustado","\n","Índices de Condición y Proporciones deVarianza","\n")
}
resul}datos %>% cor(method="pearson") %>% round(digits=2) -> corcorVerifiquemos la presencia de multicolinealidad, primero utiliaremos el diagnostico de la matriz de correlación.
corrplot(corcor, type="upper", tl.col="black")
cor(datos)## Y X1 X2 X3 X4 X5
## Y 1.0000000 0.4939233 0.03723050 0.4988930 0.42648615 0.26029349
## X1 0.4939233 1.0000000 0.20636081 0.1982719 0.37908182 0.29406162
## X2 0.0372305 0.2063608 1.00000000 -0.2476300 0.03743146 -0.08511824
## X3 0.4988930 0.1982719 -0.24762997 1.0000000 0.35967155 0.10258439
## X4 0.4264861 0.3790818 0.03743146 0.3596716 1.00000000 0.07684465
## X5 0.2602935 0.2940616 -0.08511824 0.1025844 0.07684465 1.00000000No tenemos ningun \(|\rho_{ij} |\geq 0.5\) por lo tanto no tenemos ningun indicio de multicoliealidad sin embargo sospechamos por algunos valores cercanos a este intervalo.
Utilicemos ahora el diagnostico de factores de inflación de varianza.
miscoeficientes(modelo,datos)## Coeficientes estimados, sus I.C, Vifs y Coeficientes estimados estandarizados## Estimación Limites.2.5.. Limites.97.5.. Vif Coef.Std
## (Intercept) -0.5986008691 -3.632021689 2.434819951 0.000000 0.00000000
## X1 0.2106683091 0.053437116 0.367899502 1.375666 0.30449605
## X2 0.0197512078 -0.035697988 0.075200403 1.176970 0.07487653
## X3 0.0470924905 0.020501581 0.073683400 1.270949 0.38685413
## X4 0.0105603543 -0.004080114 0.025200822 1.302309 0.15949429
## X5 0.0008995893 -0.000577038 0.002376217 1.124675 0.12518477Ya que ningun factor de inflación de varianza \(VIF_{j}> 5\) o \(VIF_{j}> 10\) no se han encontrado problemas de multicolinealidad.
Utilicemos ahora el diagnostico de índices de condición.
misDiagnostcolin(modelo)## Diagnósticos Multicolinealidad - Intercepto incluído
## Índices de Condición y Proporciones deVarianza## Val.propio cond.index Pi.intercept Pi.X1 Pi.X2 Pi.X3
## 1 5.40182563 1.000000 0.0002386136 0.0008820943 0.0002496554 0.006424605
## 2 0.30259372 4.225132 0.0001855077 0.0001147921 0.0001862061 0.155728657
## 3 0.23654019 4.778789 0.0023903752 0.0049986896 0.0038354112 0.640772446
## 4 0.03433520 12.542973 0.0223457381 0.0037191274 0.0285460562 0.114363084
## 5 0.02078651 16.120536 0.0399258551 0.9761466034 0.0293340090 0.004189488
## 6 0.00391875 37.127569 0.9349139102 0.0141386933 0.9378486621 0.078521720
## Pi.X4 Pi.X5
## 1 0.001403204 0.007901015
## 2 0.003353482 0.776605732
## 3 0.004296292 0.124324670
## 4 0.920850823 0.002073082
## 5 0.065410293 0.067694494
## 6 0.004685906 0.021401005De nuestra tabla sin centrar ya que hay valores \(10 < \sqrt{k}\leq 31.62\) y \(\sqrt{k}> 31.62\) podemos ver que en el modelo hay dos problemas de multicolinealidad moderada y uno de multicolinealidad severa.
Por ultimo utilizemos nuestro diagnostico de descomposición de varianza, en la tabla centrada.
misDiagnostcolin(modelo,centrar=T)## Diagnósticos Multicolinealidad -Intercepto ajustado
## Índices de Condición y Proporciones deVarianza## Val.propio cond.index Pi.X1 Pi.X2 Pi.X3 Pi.X4
## 1 1.7323654 1.000000 0.132673811 0.0001653996 0.10741655 0.142120663
## 2 1.2396841 1.182127 0.084395758 0.4433688520 0.13165390 0.001072787
## 3 0.9874882 1.324506 0.009577422 0.0282670390 0.05662206 0.136118366
## 4 0.5461281 1.781035 0.013604334 0.2065052224 0.64568619 0.572581500
## 5 0.4943342 1.872015 0.759748675 0.3216934870 0.05862130 0.148106685
## Pi.X5
## 1 6.484536e-02
## 2 2.244896e-06
## 3 6.375120e-01
## 4 2.662746e-03
## 5 2.949776e-01Analizando los valores \(\pi_{ij}\) de la tabla centrada podemos observar un problema de multicolinealidad moderada causado por X3 y X4 ya que hay 2 valores \(\pi_{ij}> 0.5\) asociados al mismo valor propio, de las demás relaciones no se tienen certeza de las varibles que las causan.