library(readxl)
Tabla_B4 <- read_excel("Tabla_B4.xlsx",
col_types = c("numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric"))
datos<- Tabla_B4
y<-datos$Precio_de_venta_de_casa
x1<-datos$Impuesto
x2<-datos$Cantidad_baños
x3<-datos$Tamaño_terreno
x4<-datos$Superficie_construida
x5<-datos$Cantidad_Cajones
x6<-datos$Cantidad_habitaciones
x7<-datos$Cantidad_recamaras
x8<-datos$Edad_casa
x9<-datos$Cantidad_chimeneas
datos<-data.frame(y,x1,x2,x3,x4,x5,x6,x7,x8,x9)
datos## y x1 x2 x3 x4 x5 x6 x7 x8 x9
## 1 25.9 4.9176 1.0 3.7420 0.998 1.0 7 4 42 0
## 2 29.5 5.0208 1.0 3.5310 1.500 2.0 7 4 62 0
## 3 27.9 4.5429 1.0 2.2750 1.175 1.0 6 3 40 0
## 4 25.9 4.5573 1.0 4.0500 1.232 1.0 6 3 54 0
## 5 29.9 5.0597 1.0 4.4550 1.121 1.0 6 3 42 0
## 6 29.9 3.8910 1.0 4.4550 0.988 1.0 6 3 56 0
## 7 30.9 5.8980 1.0 5.8500 1.240 1.0 7 3 51 1
## 8 28.9 5.6039 1.0 9.5200 1.501 0.0 6 3 32 0
## 9 35.9 5.8282 1.0 6.4350 1.225 2.0 6 3 32 0
## 10 31.5 5.3003 1.0 4.9883 1.552 1.0 6 3 30 0
## 11 31.0 6.2712 1.0 5.5200 0.975 1.0 5 2 30 0
## 12 30.9 5.9592 1.0 6.6660 1.121 2.0 6 3 32 0
## 13 30.0 5.0500 1.0 5.0000 1.020 0.0 5 2 46 1
## 14 36.9 8.2464 1.5 5.1500 1.664 2.0 8 4 50 0
## 15 41.9 6.6969 1.5 6.9020 1.488 1.5 7 3 22 1
## 16 40.5 7.7841 1.5 7.1020 1.376 1.0 6 3 17 0
## 17 43.9 9.0384 1.0 7.8000 1.500 1.5 7 3 23 0
## 18 37.5 5.9894 1.0 5.5200 1.256 2.0 6 3 40 1
## 19 37.9 7.5422 1.5 5.0000 1.690 1.0 7 3 22 0
## 20 44.5 8.7951 1.5 9.8900 1.820 2.0 6 4 50 1
## 21 37.9 6.0831 1.5 6.7265 1.652 1.0 6 3 44 0
## 22 38.9 8.3607 1.5 9.9150 1.777 2.0 8 4 48 1
## 23 36.9 8.1400 1.0 8.0000 1.504 2.0 7 3 3 0
## 24 45.8 9.1416 1.5 7.3262 1.831 1.5 8 4 31 0y <- datos[,c('y')]
X <- cbind(1, datos[,c('x1', 'x2','x3','x4','x5','x6','x7','x8','x9')])
Con los siguientes datos pertenecientes a el avaluo de una vivienda
Hallar la recta de regresion
Estimar \(\hat{\beta}_i\)
Estimar \(\hat{\sigma}^{2}\)
Estimar las varianza de \(\hat{\beta}_i\)
intervalo de confianza para \(\hat{\sigma}^{2}\)
Intervalos de confianza y pruebas de hipótesis para las componentes de \(\hat{\beta}\)
Usar pruebas t para evaluar la contribución de cada regresor al modelo .
**Estimaremos : \(\hat{\beta}_i\) y \(\hat{\sigma}^{2}\)
los datos de la tabla “Tabla_B4”, transformaremos a una matriz.
y<-as.matrix(y) #Matriz y
X <- as.matrix(X) #Matriz Xia).Hallaremos la traspuesta de la matriz “X”, para multiplicalo por la matriz “X” \(X^{t}X\)
Matriz_tX_X<-t(X) %*% X
round(Matriz_tX_X,4)## 1 x1 x2 x3 x4 x5 x6
## 1 24.000 153.7180 28.0000 145.8190 33.2060 31.5000 155.0000
## x1 153.718 1042.1141 185.0430 984.3218 220.0626 211.8424 1009.7381
## x2 28.000 185.0430 34.0000 174.8249 39.8550 37.5000 183.0000
## x3 145.819 984.3218 174.8249 978.7412 209.0846 197.4201 950.0014
## x4 33.206 220.0626 39.8550 209.0846 47.6987 44.9615 217.7050
## x5 31.500 211.8424 37.5000 197.4201 44.9615 49.7500 209.0000
## x6 155.000 1009.7381 183.0000 950.0014 217.7050 209.0000 1017.0000
## x7 76.000 494.3150 90.0000 466.4912 107.2130 104.0000 499.0000
## x8 899.000 5534.7107 1041.0000 5247.6232 1231.4330 1176.0000 5801.0000
## x9 6.000 40.7901 7.5000 43.0770 8.6010 8.5000 39.0000
## x7 x8 x9
## 1 76.0000 899.000 6.0000
## x1 494.3150 5534.711 40.7901
## x2 90.0000 1041.000 7.5000
## x3 466.4912 5247.623 43.0770
## x4 107.2130 1231.433 8.6010
## x5 104.0000 1176.000 8.5000
## x6 499.0000 5801.000 39.0000
## x7 248.0000 2904.000 19.0000
## x8 2904.0000 38209.000 257.0000
## x9 19.0000 257.000 6.0000b). Multiplicaremos la matriz Traspuesta de X por la matriz “y” \(X^{t}y\)
Matriz_tX_y<-t(X) %*% y
round(Matriz_tX_y,4)## [,1]
## 1 830.700
## x1 5511.926
## x2 992.850
## x3 5224.252
## x4 1176.434
## x5 1128.600
## x6 5413.300
## x7 2652.600
## x8 30344.000
## x9 223.700c). Hallaremos la matriz de los \(\hat{\beta}_i\)
\(\hat{\beta}_i\) \(=\left(X^{\ t}X\right)^{-1}\) \(\ (X^{t}\ y)\)
βi<- solve( Matriz_tX_X)%*% (Matriz_tX_y)
round(βi,4)## [,1]
## 1 17.5071
## x1 2.1146
## x2 6.6573
## x3 -0.0778
## x4 3.3929
## x5 1.8655
## x6 -1.0829
## x7 -0.7522
## x8 -0.0519
## x9 1.7737Entonces la ecuacion de regresion sera:
B<-round(βi,4)
paste( "βi= ",B[1],"+" ,B[2],"Impuesto +",B[3],"Cantidad de baños " ,B[4],"Iamaño del terreno +" ,B[5],"Superficie construida +",B[6],"Cantidad de cajones ",B[7],"Cantidad de habitaciones ",B[8],"Cantidad de recamaras",B[9],"Edad de la casa +",B[10],"Cantidad de chimeneas")## [1] "ßi= 17.5071 + 2.1146 Impuesto + 6.6573 Cantidad de baños -0.0778 Iamaño del terreno + 3.3929 Superficie construida + 1.8655 Cantidad de cajones -1.0829 Cantidad de habitaciones -0.7522 Cantidad de recamaras -0.0519 Edad de la casa + 1.7737 Cantidad de chimeneas"Ahora estimaremos \(\hat{\sigma}^{2}\)
Usaremos la siguiente formula:
\(\hat{\sigma}^{2}=\frac{1}{n-(p+1)}(\mathbf{y}-\mathbf{X} \hat{\beta})^{\prime}(\mathbf{y}-\mathbf{X} \hat{\beta})\)
s2h <- t(y - X%*%βi) %*% (y - X%*%βi) / (24 - 10)
s2h <- as.numeric(s2h); round(s2h, 4)## [1] 8.3912Dando como resultado: \(\hat{\sigma}^{2}=\frac{1}{n-(p+1)}(\mathbf{y}-\mathbf{X} \hat{\beta})^{\prime}(\mathbf{y}-\mathbf{X} \hat{\beta})=8.3912\)
Ahora estimaremos la varianza de \(\hat{\beta}_i\)
Varianza.B <- s2h * solve(t(X) %*% X); round(Varianza.B, 4)## 1 x1 x2 x3 x4 x5 x6 x7 x8
## 1 34.2193 0.2689 -4.8747 -1.0579 2.8253 0.2457 -5.3608 4.1483 -0.1855
## x1 0.2689 0.7753 -1.3094 -0.2102 -0.7407 -0.4459 -0.4098 0.3572 0.0182
## x2 -4.8747 -1.3094 16.5075 0.4979 -7.4824 1.3380 0.6993 -1.4136 -0.0018
## x3 -1.0579 -0.2102 0.4979 0.2560 -0.6385 0.1307 0.2634 -0.2877 0.0075
## x4 2.8253 -0.7407 -7.4824 -0.6385 18.4310 0.3872 -0.8428 -1.8207 -0.0141
## x5 0.2457 -0.4459 1.3380 0.1307 0.3872 1.8436 0.2064 -1.3527 0.0098
## x6 -5.3608 -0.4098 0.6993 0.2634 -0.8428 0.2064 1.8873 -2.0365 0.0226
## x7 4.1483 0.3572 -1.4136 -0.2877 -1.8207 -1.3527 -2.0365 5.4519 -0.0858
## x8 -0.1855 0.0182 -0.0018 0.0075 -0.0141 0.0098 0.0226 -0.0858 0.0046
## x9 2.8960 0.0515 -1.5938 -0.3659 1.0273 -0.4709 -0.4899 1.4302 -0.0548
## x9
## 1 2.8960
## x1 0.0515
## x2 -1.5938
## x3 -0.3659
## x4 1.0273
## x5 -0.4709
## x6 -0.4899
## x7 1.4302
## x8 -0.0548
## x9 2.9704Ahora hallaremos el intevalo de confianza para \(\hat{\sigma}^{2}\)
Un intervalo de confianza \(100\ (1-\alpha)\%\) para \(\hat{\sigma}^{2}\) esta dado por:
\(\left(\frac{(n-p-1) \sigma_{M C O}^{2}}{\chi_{n-p-1}^{2}(\alpha / 2)}, \frac{(n-p-1) \sigma_{M C O}^{2}}{\chi_{n-p-1}^{2}(1-\alpha / 2)}\right)\)
Calculamos el intervalo de confianza \(95%\) para \(\hat{\sigma}^{2}\)
n<-24
gl<-10
gamma1 <- qchisq(0.025, n - gl)
gamma2 <- qchisq(0.975, n - gl)
s2_LI <- (n - gl) * s2h / gamma2; s2_LI## [1] 4.497749s2_LS <- (n - gl) * s2h / gamma1; s2_LS## [1] 20.87088paste("El limite inferior con 95% de confianza es: ",round(s2_LI,4))## [1] "El limite inferior con 95% de confianza es: 4.4977"paste("El limite superior con 95% de confianza es: ",round(s2_LS,4))## [1] "El limite superior con 95% de confianza es: 20.8709"Intervalos de confianza y pruebas de hipótesis para las componentes de \(\hat{\beta}\)
Los intervalos de confianza para los componentes de \(\hat{\beta}\) están dados por:
\(\hat{\beta}{i} \pm t_{n-p-1}(\alpha / 2) \sqrt{\hat{V}(\hat{\boldsymbol{\beta}})_{i i}}\)
donde \(\hat{\beta}_i\) es la \(i-esima\) entrada de \(\hat{\beta}\), \(t_{n-p-1}(\alpha / 2)\) es el cuantil superior \((\alpha / 2)\) de una distribucion \(t_{n-p-1}\) y \(\hat{V}(\hat{\boldsymbol{\beta}})_{i i}\) es el \(i-esimo\) elemento de la diagonal de \(\hat{V}(\hat{\boldsymbol{\beta}})_{i i}\).
Para contrastar las hipotesis
\(H_0: \beta_{i}=0\)
\(H_1: \beta_{i}\neq 0\)
Para ello se utiliza el estadistico
\(t_{i}=\frac{\hat{\beta}{i}-b}{\sqrt{\hat{V}(\hat{\beta}){i i}}}\)
la regla de decision es rechazar si \(\left|t_{i}\right|>t_{n-p-1}(\alpha / 2)\)
Para \(\beta_0\)
b0_LI <- B[1] - qt(0.975, n - gl) * sqrt(Varianza.B[1,1]);
b0_LS <- B[1] + qt(0.975, n - gl) * sqrt(Varianza.B[1,1]);
paste("El limite inferior es: ",b0_LI,
"El limite superior es: ",b0_LS)## [1] "El limite inferior es: 4.9606799590982 El limite superior es: 30.0535200409018"t0 <- abs(B[1]/sqrt(Varianza.B[1,1])); t0 # Estadístico## [1] 2.992806tt <- qt(0.975, n - gl); tt # Cuantil de la distribución t ## [1] 2.144787Conclusion:
Como |t1|>tn−p−1(0.025), se rechaza la hipótesis nula con una significancia α=0.05, es decir, la evidencia sostiene que el modelo incluya a β0.
Para \(\beta_1\)
b0_LI <- B[2] - qt(0.975, n - gl) * sqrt(Varianza.B[2,2]);
b0_LS <- B[2] + qt(0.975, n - gl) * sqrt(Varianza.B[2,2]);
paste("El limite inferior es: ",b0_LI,
"El limite superior es: ",b0_LS)## [1] "El limite inferior es: 0.226045518628881 El limite superior es: 4.00315448137112"t0 <- abs(B[2]/sqrt(Varianza.B[2,2])); t0 # Estadístico## [1] 2.401501tt <- qt(0.975, n - gl); tt # Cuantil de la distribución t ## [1] 2.144787Conclusion:
Como |t1|>tn−p−1(0.025), se rechaza la hipótesis nula con una significancia α=0.05, es decir, la evidencia sostiene que el modelo incluya a β1.
Para \(\beta_2\)
b0_LI <- B[3] - qt(0.975, n - gl) * sqrt(Varianza.B[3,3]);
b0_LS <- B[3] + qt(0.975, n - gl) * sqrt(Varianza.B[3,3]);
paste("El limite inferior es: ",b0_LI,
"El limite superior es: ",b0_LS)## [1] "El limite inferior es: -2.05685645246859 El limite superior es: 15.3714564524686"t0 <- abs(B[3]/sqrt(Varianza.B[3,3])); t0 # Estadístico## [1] 1.638539tt <- qt(0.975, n - gl); tt # Cuantil de la distribución t ## [1] 2.144787Conclusion:
Como |t1|<tn−p−1(0.025), no se rechaza la hipótesis nula con una significancia α=0.05, es decir, la evidencia sostiene que el modelo no debería incluir a β2.
Para \(\beta_3\)
b0_LI <- B[4] - qt(0.975, n - gl) * sqrt(Varianza.B[4,4]);
b0_LS <- B[4] + qt(0.975, n - gl) * sqrt(Varianza.B[4,4]);
paste("El limite inferior es: ",b0_LI,
"El limite superior es: ",b0_LS)## [1] "El limite inferior es: -1.16300113500484 El limite superior es: 1.00740113500484"t0 <- abs(B[4]/sqrt(Varianza.B[4,4])); t0 # Estadístico## [1] 0.1537636tt <- qt(0.975, n - gl); tt # Cuantil de la distribución t ## [1] 2.144787Conclusion:
Como |t1|<tn−p−1(0.025), no se rechaza la hipótesis nula con una significancia α=0.05, es decir, la evidencia sostiene que el modelo no debería incluir a β3.
Para \(\beta_4\)
b0_LI <- B[5] - qt(0.975, n - gl) * sqrt(Varianza.B[5,5]);
b0_LS <- B[5] + qt(0.975, n - gl) * sqrt(Varianza.B[5,5]);
paste("El limite inferior es: ",b0_LI,
"El limite superior es: ",b0_LS)## [1] "El limite inferior es: -5.81494631934622 El limite superior es: 12.6007463193462"t0 <- abs(B[5]/sqrt(Varianza.B[5,5])); t0 # Estadístico## [1] 0.7903093tt <- qt(0.975, n - gl); tt # Cuantil de la distribución t ## [1] 2.144787Conclusion:
Como |t1|<tn−p−1(0.025), no se rechaza la hipótesis nula con una significancia α=0.05, es decir, la evidencia sostiene que el modelo no debería incluir a β4.
Para \(\beta_5\)
b0_LI <- B[6] - qt(0.975, n - gl) * sqrt(Varianza.B[6,6]);
b0_LS <- B[6] + qt(0.975, n - gl) * sqrt(Varianza.B[6,6]);
paste("El limite inferior es: ",b0_LI,
"El limite superior es: ",b0_LS)## [1] "El limite inferior es: -1.0466467213881 El limite superior es: 4.7776467213881"t0 <- abs(B[6]/sqrt(Varianza.B[6,6])); t0 # Estadístico## [1] 1.373935tt <- qt(0.975, n - gl); tt # Cuantil de la distribución t ## [1] 2.144787Conclusion:
Como |t1|<tn−p−1(0.025), no se rechaza la hipótesis nula con una significancia α=0.05, es decir, la evidencia sostiene que el modelo no debería incluir a β5.
Para \(\beta_6\)
b0_LI <- B[7] - qt(0.975, n - gl) * sqrt(Varianza.B[7,7]);
b0_LS <- B[7] + qt(0.975, n - gl) * sqrt(Varianza.B[7,7]);
paste("El limite inferior es: ",b0_LI,
"El limite superior es: ",b0_LS)## [1] "El limite inferior es: -4.02935177577978 El limite superior es: 1.86355177577978"t0 <- abs(B[7]/sqrt(Varianza.B[7,7])); t0 # Estadístico## [1] 0.7882666tt <- qt(0.975, n - gl); tt # Cuantil de la distribución t ## [1] 2.144787Conclusion:
Como |t1|<tn−p−1(0.025), no se rechaza la hipótesis nula con una significancia α=0.05, es decir, la evidencia sostiene que el modelo no debería incluir a β6.
Para \(\beta_7\)
b0_LI <- B[8] - qt(0.975, n - gl) * sqrt(Varianza.B[8,8]);
b0_LS <- B[8] + qt(0.975, n - gl) * sqrt(Varianza.B[8,8]);
paste("El limite inferior es: ",b0_LI,
"El limite superior es: ",b0_LS)## [1] "El limite inferior es: -5.7601124448246 El limite superior es: 4.2557124448246"t0 <- abs(B[8]/sqrt(Varianza.B[8,8])); t0 # Estadístico## [1] 0.3221519tt <- qt(0.975, n - gl); tt # Cuantil de la distribución t ## [1] 2.144787Conclusion:
Como |t1|<tn−p−1(0.025), no se rechaza la hipótesis nula con una significancia α=0.05, es decir, la evidencia sostiene que el modelo no debería incluir a β7.
Para \(\beta_8\)
b0_LI <- B[9] - qt(0.975, n - gl) * sqrt(Varianza.B[9,9]);
b0_LS <- B[9] + qt(0.975, n - gl) * sqrt(Varianza.B[9,9]);
paste("El limite inferior es: ",b0_LI,
"El limite superior es: ",b0_LS)## [1] "El limite inferior es: -0.196857908307772 El limite superior es: 0.0930579083077721"t0 <- abs(B[9]/sqrt(Varianza.B[9,9])); t0 # Estadístico## [1] 0.7679086tt <- qt(0.975, n - gl); tt # Cuantil de la distribución t ## [1] 2.144787Conclusion:
Como |t1|<tn−p−1(0.025), no se rechaza la hipótesis nula con una significancia α=0.05, es decir, la evidencia sostiene que el modelo no debería incluir a β8.
Para \(\beta_9\)
b0_LI <- B[10] - qt(0.975, n - gl) * sqrt(Varianza.B[10,10]);
b0_LS <- B[10] + qt(0.975, n - gl) * sqrt(Varianza.B[10,10]);
paste("El limite inferior es: ",b0_LI,
"El limite superior es: ",b0_LS)## [1] "El limite inferior es: -1.92282368680731 El limite superior es: 5.47022368680731"t0 <- abs(B[10]/sqrt(Varianza.B[10,10])); t0 # Estadístico## [1] 1.029131tt <- qt(0.975, n - gl); tt # Cuantil de la distribución t ## [1] 2.144787Conclusion:
Como |t1|<tn−p−1(0.025), no se rechaza la hipótesis nula con una significancia α=0.05, es decir, la evidencia sostiene que el modelo no debería incluir a β9.
Analisis
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericregresion.lineal.multiple<-lm(y~x1+x2+x3+x4+x5+x6+x7+x8+x9)
summary(regresion.lineal.multiple)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8730 -1.8497 -0.3924 1.7597 4.3724
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.50709 5.84973 2.993 0.00969 **
## x1 2.11456 0.88053 2.401 0.03078 *
## x2 6.65733 4.06295 1.639 0.12358
## x3 -0.07781 0.50597 -0.154 0.87998
## x4 3.39286 4.29313 0.790 0.44253
## x5 1.86553 1.35778 1.374 0.19106
## x6 -1.08286 1.37377 -0.788 0.44370
## x7 -0.75222 2.33492 -0.322 0.75209
## x8 -0.05187 0.06759 -0.768 0.45553
## x9 1.77366 1.72349 1.029 0.32088
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.897 on 14 degrees of freedom
## Multiple R-squared: 0.8583, Adjusted R-squared: 0.7672
## F-statistic: 9.422 on 9 and 14 DF, p-value: 0.0001467Hallamos el \(R^2\) y \(R^2_{Adj}\)
\(R^2\)
R2<-summary(regresion.lineal.multiple)$r.squared
paste("El valor del coeficiente de detreminacion es:",round(R2*100,2),"%" )## [1] "El valor del coeficiente de detreminacion es: 85.83 %"\(R^2_{Adj}\)
R2.ajustado<-summary(regresion.lineal.multiple)$adj.r.squared
paste("El valor del coeficiente de detreminacion ajustado es:",round(R2.ajustado*100,2),"%" )## [1] "El valor del coeficiente de detreminacion ajustado es: 76.72 %"