Lectura de la base de datos

require(ggplot2)
## Loading required package: ggplot2
datos2 = read.csv("~/Descargas/datos2.txt", sep="")
source("~/Descargas/functions.R")

Modelo de regresión lineal

mod = lm(Y~., data=datos2)
summary(mod)
## 
## Call:
## lm(formula = Y ~ ., data = datos2)
## 
## 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-06

Para tener una mejor apreciación del modelo se presenta estadísticos adicionales

anova(mod)
## Analysis of Variance Table
## 
## Response: Y
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## X1         1 27.872 27.8721 26.0185 3.770e-06 ***
## X2         1  0.499  0.4995  0.4662    0.4974    
## X3         1 19.028 19.0275 17.7621 8.703e-05 ***
## X4         1  2.054  2.0542  1.9176    0.1713    
## X5         1  1.592  1.5919  1.4861    0.2277    
## Residuals 59 63.203  1.0712                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Significancia de la regresión

Al analizar los resultados obtenidos por el resumen del modelo, se puede notar que el valor-p de la regresión está muy cerca a cero, con lo cual podemos concluir que la regresión sí es significativa.

Significación de los parámetros

Similarmente se utilizará el valor-p para determinar la significancia de estos. Así, se puede notar que para los parámetros \(X1\) y \(X3\) estos valores son casi cero y, por ende, se dice que sí afectan al modelo. En cambio, los parámetros \(X2\), \(X4\) y \(X5\) tiene un valor-p suficientemente grande y, entonces, se concluye que no son significantes para el modelo.

Intrepretación de los parámetros

  • En \(X1\) se puede decir que por cada día que pasa el paciente en el hospital, la probabilidad de adquirir la infección aumenta en aproximadamente un \(21\%\).
  • A medida en que aumenta la razón de los cultivos realizados en pacientes sin síntomas, el riesgo de infección aumenta en aproximadamente \(1\%\)
  • Si se incrementa en uno el número de camas en el hospital, la probabilidad de infectarse aumentará en un \(4%\) aproximadamente.
  • Si la cantidad de pacientes promedio en el hospital incrementa, entonces el riesgo de adquirir la infección aumenta en \(1\%\) aproximadamente.
  • Finalmente, si se incrementa la cantidad de enfermeras en el hospital, la probabilidad de adquirir la infección no incrementa significativamente.

Interpretación del \(R^2\)

Del resumen estadístico del modelo, se puede extraer el \(R^2\). En este caso se tiene \(R^2 = 0.4468\), lo que explica que el modelo es capaz de predecir de forma adecuada el \(44\%\) de la veces la probabilidad de contraer la infección.

Significancia simultánea

Con el fin de analizar cuán importante son las variables asociadas a los parámetros \(\beta_i\) con mayor valor-p (\(X_2, X_4\,y \,X_5\)), se procede a realizar la regresión lineal en cual se involucran estas variables.

# myAllRegTable(mod)

Para saber la significancia de este modelo se plantea la siguiente prueba de hipótesis: \(H_0: X_2 = 0, X_4 = 0, X_5 = 0\) vs \(H_1: X_i \ne 0, i = 2,4,5\)

SSE_mod = 63.203
MSE_mod = 1.0712
SSE_modred = 87.321
SSR_extra = SSE_modred - SSE_mod

F0 = (SSR_extra/3) / MSE_mod
F_critico = qf(0.95, df1 = 3, df2 = 59)
F0 > F_critico
## [1] TRUE

Debido a que \(F_0 > F\_critico\) se puede rechazar la hipótesis nula (\(H_0\)), es decir, alguna de las variables sí es significativa para el modelo y, por tanto, no se pueden descartar.

Hipótesis Lineal

Se quiere probar lo siguiente: \(H_0 : X_1 = X_4, X_2 = 0\) vs \(H_1: X_1 \ne X_4\,ó\,X_2 \ne 0\). Si se reescribe lo anterior, se tiene: \(H_0: X_1 - X_4 = 0, X_2 = 0\). Esto expresado de forma matricial queda así: \[H_0: \begin{bmatrix} 0 & 1 & 0 & 0 & -1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \\ \beta_4 \\ \beta_5 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}\] donde \(\textbf L = \begin{bmatrix} 0 & 1 & 0 & 0 & -1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix}\).
Así, el modelo reducido queda expresado de la siguiente forma: \(Y = \beta_0 + \beta_1(X_1 + X_4) + \beta_3X_3 + \beta_5X_5\). Finalmente, el estadístico de prueba se define así:

SSE_mod = 63.203
MSE_mod = 1.0712
SSE_modred = 87.321
SSH = SSE_modred - SSE_mod
MSH = SSH / 2
F0 = MSH / MSE_mod

Validación de supuestos

Media cero

Para validar este supuesto, se obtienen los residuales del modelo y queda claro que la media de estos es cero, por tanto se cumple por defecto.

ei=mod$residuals

round(mean(ei),3)
## [1] 0

Varianza constante

Al graficar los residuales vs los ajustados se puede notar que no hay una tendencia fuerte que indique problemas. Se valida por medio de la gráfica.

yi_mod=mod$fitted.values
plot(yi_mod,ei,main="Varianza Constante (ei)")

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'

Normalidad

En esta se plantea una prueba de Shapiro-Wilk como sigue: \(H_0: Normalidad\) vs \(H_1: No\,\,normalidad\). De la cual se observa que \(p\_value = 0.5979\), por tanto a un nivel de significancia de \(\alpha=0.05\) se puede concluir que no se puede rechazar la hipótesis nula \(H_0\) ya que \(p\_value \nless alpha\), es decir, los errores sí tienen una distribución normal.

hist(ei,col="gray")

qqnorm(ei)
qqline(ei,col="red")

shapiro.test(ei)
## 
##  Shapiro-Wilk normality test
## 
## data:  ei
## W = 0.98462, p-value = 0.5979

Independencia

Ya que los registros no están asociados a datos en el tiempo, entonces no se tiene un orden temporal entre estos. Se valida por la definición de los datos.

Multicolinealidad

Para detectar posible relaciones fuertes entre las variables, se presenta la matrix de Pearson con el fin de relacionar todas contra todas.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
datos2 %>% glimpse()
## Rows: 65
## Columns: 6
## $ Y  <dbl> 3.7, 2.8, 4.2, 6.2, 5.7, 4.5, 1.6, 5.1, 4.1, 4.4, 5.0, 4.3, 5.3, 4.…
## $ X1 <dbl> 7.58, 9.97, 8.88, 10.15, 11.18, 9.61, 8.82, 10.30, 10.47, 10.02, 9.…
## $ X2 <dbl> 56.7, 58.2, 51.5, 51.9, 51.0, 52.4, 58.2, 59.6, 53.2, 49.5, 52.3, 4…
## $ X3 <dbl> 20.8, 16.5, 10.1, 16.4, 18.8, 6.9, 3.8, 27.8, 5.7, 8.3, 17.6, 16.4,…
## $ X4 <dbl> 88.0, 76.5, 86.9, 59.2, 55.9, 87.2, 51.7, 88.9, 69.1, 93.0, 95.9, 6…
## $ X5 <int> 97, 90, 305, 568, 595, 487, 80, 175, 196, 265, 270, 318, 99, 600, 2…
datos2 %>% cor(method="pearson") %>% round(digits=2) -> mat_cor
mat_cor
##       Y   X1    X2    X3   X4    X5
## Y  1.00 0.49  0.04  0.50 0.43  0.26
## X1 0.49 1.00  0.21  0.20 0.38  0.29
## X2 0.04 0.21  1.00 -0.25 0.04 -0.09
## X3 0.50 0.20 -0.25  1.00 0.36  0.10
## X4 0.43 0.38  0.04  0.36 1.00  0.08
## X5 0.26 0.29 -0.09  0.10 0.08  1.00
library(corrplot)
## corrplot 0.90 loaded
corrplot(mat_cor, type="upper", order="hclust", tl.col="black", tl.srt=45)

De lo anterior se puede notar que no existe una relación entre las variables predecitoras ya que no existe algún par que supere \(|0.5|\) en su coefieciente de correlación.