Importamos los datos desde un hoja de excel

library(readxl)
EJERCICIO_3_14_1_ <- read_excel("EJERCICIO 3.14(1).xlsx", 
    sheet = "Datos para R", col_types = c("numeric", 
        "numeric", "numeric"))
View(EJERCICIO_3_14_1_)
datos<-EJERCICIO_3_14_1_

1. Ajustar un modelo de regresión lineal múltiple que relacione la viscosidad con los dos regresores

Converitmos en matriz los datos

y <- datos[,c('y')] # Convertimos los dato de "y" en una matriz.
X <- cbind(1, datos[,c('x1', 'x2')])  # Convertimos los datos de "X" en una matriz.
y<-as.matrix(y) # Matriz y.
X <- as.matrix(X) # Matriz Xi.

a).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
## 1    40.000  26.0310  1400.000
## x1   26.031  18.6746   911.085
## x2 1400.000 911.0850 82000.000

b). Multiplicaremos la matriz Traspuesta de X por la matriz “y” \(X^{t}y\) 

Matriz_tX_y<-t(X) %*% y
round(Matriz_tX_y,4)
##           y
## 1   41.9314
## x1  29.7286
## x2 951.8470

c). 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)
##          y
## 1   0.6794
## x1  1.4073
## x2 -0.0156

Entonces la ecuacion de regresion sera:

B<-round(βi,4)
paste( "Viscosidad cinematica= ",B[1],"+" ,B[2],"(relacion de 2-metoxietanos a 1,2-dimetoxietano) ",B[3]," (Temperatura) ")
## [1] "Viscosidad cinematica=  0.6794 + 1.4073 (relacion de 2-metoxietanos a 1,2-dimetoxietano)  -0.0156  (Temperatura) "

2.Probar la significancia de la regresión. ¿Qué conclusiones se pueden sacar?

Para ello se debe probr que:

\(H_{0}=\beta_{0}=\beta_{1}=\beta_{2}=0\)

\(H_{1} \neq \beta_{j}\)

Primero hallaremos la suma total de cuadrado: \((SS_T)\)

\(S S_{T}=\sum y_{i}^{2}-\frac{\left(\sum y_{i}\right)^{2}}{n}\)

y.cua<-y^2
n<-40
SSt=sum(y.cua)-((sum(y)^2)/n);SSt
## [1] 13.98399

Enntonces:

\(S S_{T}=\sum y_{i}^{2}-\frac{\left(\sum y_{i}\right)^{2}}{n}= 13.98399\)

Ahora hallaremos la suma de cuadrados de la regresión: \((SS_R)\)

\(S S_{R}=\hat{\beta}^{\prime} X^{\prime} y-\frac{\left(\sum y_{i}\right)^{2}}{n}\)

SSr<-(t(B)%*%Matriz_tX_y)-((sum(y)^2)/n);SSr
##          y
## y 11.52044

Entonces:

\(S S_{R}=\hat{\beta}^{\prime} X^{\prime} y-\frac{\left(\sum y_{i}\right)^{2}}{n}=11.52045\)

Ahora hallaremos la suma de cuadrados residuales \((SS_{Res}=SS_T-SS_R=2.46355)\) 

SSres<-SSt-SSr;SSres
##          y
## y 2.463544

Para probar que \(H_{0}:\beta_{0}=\beta_{1}=\beta_{2}=0\), se calcula el estadistico siguiente:

\(F_{0}=\frac{\frac{S S_{R}}{k}}{\frac{S S_{R e s}}{n-k-1}}\)

k<-2
Fo<-(SSr/k)/(SSres/(n-k-1));Fo
##          y
## y 86.51282

Hallamos el F tabulado

F.tab<-qf(0.05,k,37,lower.tail = F);F.tab
## [1] 3.251924

Desicion

ifelse(Fo>F.tab,"Se rechaza la hipotesis nula(la regresion es significativa)","Se acepta la hipotesis nula")
##   y                                                            
## y "Se rechaza la hipotesis nula(la regresion es significativa)"

3. Hacer pruebas t para evaluar la contribución de cada regresor al modelo. Comentar los resultados

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})\)

gl<-3
s2h <- t(y - X%*%βi) %*% (y - X%*%βi) / (n - gl)
s2h <- as.numeric(s2h); round(s2h, 4)
## [1] 0.0673

Dando como resultado: \(\hat{\sigma}^{2}=\frac{1}{n-(p+1)}(\mathbf{y}-\mathbf{X} \hat{\beta})^{\prime}(\mathbf{y}-\mathbf{X} \hat{\beta})=0.0673\)

Ahora estimaremos la varianza de \(\hat{\beta}_i\) 

Varianza.B <- s2h * solve(t(X) %*% X); round(Varianza.B, 4)
##          1      x1     x2
## 1   0.0206 -0.0252 -1e-04
## x1 -0.0252  0.0388  0e+00
## x2 -0.0001  0.0000  0e+00

Para \(\beta_{0}\) 

gl<-3
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:  0.388576388769405 El limite superior es:  0.970223611230595"
t0 <- abs(B[1]/sqrt(Varianza.B[1,1])); t0 # Estadístico
## [1] 4.733437
tt <- qt(0.975, n - gl); tt   # Cuantil de la distribución t 
## [1] 2.026192

Desicion

ifelse(t0>tt,"Se rechaza la hipotesis nula(la regresion es significativa)","Se acepta la hipotesis nula")
## [1] "Se rechaza la hipotesis nula(la regresion es significativa)"

Conclusion:

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:  1.00829142341049 El limite superior es:  1.80630857658951"
t0 <- abs(B[2]/sqrt(Varianza.B[2,2])); t0 # Estadístico
## [1] 7.146364
tt <- qt(0.975, n - gl); tt   # Cuantil de la distribución t 
## [1] 2.026192

Desicion

ifelse(t0>tt,"Se rechaza la hipotesis nula(la regresion es significativa)","Se acepta la hipotesis nula")
## [1] "Se rechaza la hipotesis nula(la regresion es significativa)"

Conclusion:

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:  -0.018492594510196 El limite superior es:  -0.012707405489804"
t0 <- abs(B[3]/sqrt(Varianza.B[3,3])); t0 # Estadístico
## [1] 10.92742
tt <- qt(0.975, n - gl); tt   # Cuantil de la distribución t 
## [1] 2.026192

Desicion

ifelse(t0>tt,"Se rechaza la hipotesis nula(la regresion es significativa)","Se acepta la hipotesis nula")
## [1] "Se rechaza la hipotesis nula(la regresion es significativa)"

Conclusion:

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 β2.

4. Calcular \(R^2\) y \(R^2_{Adj}\) para este modelo. Comparar esos valores con los de \(R^2\) y \(R^2_{Adj}\) para el modelo de regresión lineal simple que relacione la viscosidad sólo con la temperatura. Comentar los resultados obtenidos.

Calcularemos \(R^2\) y \(R^2_{Adj}\) del modelo multilineal

\(R^2=\frac{SS_R}{SS_T}\)

R2.multi<-SSr/SSt;R2.multi
##          y
## y 0.823831

El valor de

\(R^2=82.3830681\)%

Ahora calculamos

\(R^2_{Adj}\)

\(R^2_{Adj}=1-\frac{SS_{Res}/(n-p)}{SS_T /(n-1)}\)

R2.ajus.multi<-1-((SSres/(n-3))/(SSt/(n-1)))
R2.ajus.multi
##           y
## y 0.8143084

El valor de

\(R^2_{Adj}=81.9194999\)%

Calcularemos \(R^2\) y \(R^2_{Adj}\) del modelo linael

modelo.lineal<-lm(datos$y~datos$x2)
summary(modelo.lineal)
## 
## Call:
## lm(formula = datos$y ~ datos$x2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.60279 -0.27348 -0.02503  0.19194  1.37642 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.595295   0.098404  16.212  < 2e-16 ***
## datos$x2    -0.015629   0.002173  -7.191 1.36e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3948 on 38 degrees of freedom
## Multiple R-squared:  0.5764, Adjusted R-squared:  0.5653 
## F-statistic: 51.71 on 1 and 38 DF,  p-value: 1.358e-08

El modelo de regresion es: $ y=1.595295-0.015629 x_2$

R2.line<-(cor(datos$y,datos$x2))^2;R2.line
## [1] 0.5764172
R2.line.ajus<-summary(modelo.lineal)$adj.r.squared

El valor de

\(R^2=-57.6417134\)%

El valor de

\(R^2_{Adj}=56.53\)%

se<-sqrt(s2h*(solve(Matriz_tX_X))[3,3]);se
## [1] 0.001427601
to2<-B[3]/se;to2
## [1] -10.92742

Analisis de Datos 

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
regresion.lineal.multiple<-lm(datos$y~datos$x1+datos$x2)
summary(regresion.lineal.multiple)
## 
## Call:
## lm(formula = datos$y ~ datos$x1 + datos$x2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22179 -0.18102 -0.08439  0.09111  0.99908 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.679439   0.143532   4.734 3.20e-05 ***
## datos$x1     1.407331   0.196925   7.147 1.81e-08 ***
## datos$x2    -0.015629   0.001428 -10.948 3.67e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2593 on 37 degrees of freedom
## Multiple R-squared:  0.822,  Adjusted R-squared:  0.8124 
## F-statistic: 85.46 on 2 and 37 DF,  p-value: 1.351e-14