Taller / Regresión Lineal simple

Introducción

El modelo de regresión lineal simple o ajuste lineal en el presente curso tiene dos objetivos, el primero es hallar una relación de dependencia etre dos variables que presenta el modelo X e Y, donde denominamos a la variable dependiente a Y, y las variables independientes a X. Por otra parte, el segundo objetivo es el poder predecir, es decir, dado un valor de la variable X se podrá conocer el valor que podría tener un valor para la variable Y.

\[Y = \beta_0 + \beta_1 X +\epsilon\]

Existen además hipótesis básicas del modelo estas son:

• Linealidad: consiste en que la media de variable dependiente toma un valor inicial \(\beta_0\) cuando la variable explicativa X es cero, además de que la media crece en una cantidad fija de \(\beta_1\) cada vez que se incrementa en una unidad a X.
• Homocedasticidad: nos indica que la varianza del error \(\epsilon\) es la misma cualquiera sea el valor de la variable explicativa X. \[V(\epsilon/X=x)=\sigma^2 \hspace{1em} \forall x.\]
• Normalidad: indica que el error del modelo \(\epsilon\) tiene distribución normal \[\epsilon \in N (0,\sigma^2)\]
• Independencia:esta suposición dice que los n errores de una muestra no observados serían mutuamente independientes.

En el presente documento utilizaremos funciones de R para probar las hipòtesis explicadas anteriormente mediante el siguiente ejemplo:

Ejemplo

En este tema consideraremos el ejemplo de Sheather (2009) sobre tiempos en una cadena de producción. Se dispone del tiempo (en minutos) que lleva producir cada pedido, que denotamos por Y , junto con el número de unidades solicitadas, denotado por X, para 20 pedidos seleccionados aleatoriamente. Los datos figuran a continuación:

#Adecuacion del modelo
production <- read.table("production.txt",header=TRUE)
attach(production)
production

##    Case RunTime RunSize
## 1     1     195     175
## 2     2     215     189
## 3     3     243     344
## 4     4     162      88
## 5     5     185     114
## 6     6     231     338
## 7     7     234     271
## 8     8     166     173
## 9     9     253     284
## 10   10     196     277
## 11   11     220     337
## 12   12     168      58
## 13   13     207     146
## 14   14     225     277
## 15   15     169     123
## 16   16     215     227
## 17   17     147      63
## 18   18     230     337
## 19   19     208     146
## 20   20     172      68

summary(production$RunTime)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   147.0   171.2   207.5   202.1   226.2   253.0

summary(production$RunSize)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    58.0   120.8   182.0   201.8   278.8   344.0

m1 <- lm(RunTime~RunSize)
summary(m1)

## 
## Call:
## lm(formula = RunTime ~ RunSize)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.597 -11.079   3.329   8.302  29.627 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 149.74770    8.32815   17.98 6.00e-13 ***
## RunSize       0.25924    0.03714    6.98 1.61e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.25 on 18 degrees of freedom
## Multiple R-squared:  0.7302, Adjusted R-squared:  0.7152 
## F-statistic: 48.72 on 1 and 18 DF,  p-value: 1.615e-06

m1$coefficients

## (Intercept)     RunSize 
## 149.7477030   0.2592431

anova(m1)

## Analysis of Variance Table
## 
## Response: RunTime
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## RunSize    1 12868.4 12868.4  48.717 1.615e-06 ***
## Residuals 18  4754.6   264.1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hipótesis de Linealidad

Vamos a realizar la prueba de hipótesis para la Linealidad \[H_0: Y = \beta_0 + \beta_1 X \hspace{1em} \text{Es lineal}\] \[H_a: Y = \beta_0 + \beta_1 X \hspace{1em} \text{No es lineal}\]

Y_reg<-function(x){
return(as.numeric(m1$coefficients[1])+as.numeric(m1$coefficients[2])*as.numeric(x))
}

Test de Linealidad:

Prueba_de_linealidad<-function(X,Y,regresion_lineal){
  Yt<-regresion_lineal(X)
  F_obs<- sum((Yt-mean(Y))^2)*(length(X)-2)/sum((Y-Yt)^2)
  pval <- 1-pf(F_obs,lower.tail = TRUE,df1 = 1,df2 = (length(X)-2))
  if (pval > 0.05){
    resultado<-TRUE
    h<-c("No se rechaza Ho.")
  }
  else{
    resultado<-FALSE
    h<-c("Se rechaza Ho.")
  }
  return(matrix(c("p_value","RESULTADO",pval,h),2,2))
}
Prueba_de_linealidad(RunSize,RunTime,Y_reg)

##      [,1]        [,2]                 
## [1,] "p_value"   "1.6148628918522e-06"
## [2,] "RESULTADO" "Se rechaza Ho."

Hipótesis de Homocedasticidad

\[H_0: \sigma_1=\sigma_2\] \[H_1: \sigma_1\not=\sigma_2\] Hemos agrupado los datos con el fin de contrastar las varianzas en cada uno de estos conjuntos.

Prueba_homocedasticidad<-function(X,Y,regresion_lineal)
{dist<- function(X,Y,regresion_lineal){
 yreg<- regresion_lineal(X)
 return(Y-yreg)
  }
  error<- dist(RunSize,RunTime,Y_reg)
  p<-bartlett.test(error,rep(c("Parte_1","Parte_2"),each = 10))
  pval<-p$p.value
  if (pval > 0.05){
    resultado<-TRUE
    h<-c("No se rechaza Ho")
  }
  else{
    resultado<-FALSE
    h<-c("Se rechaza Ho")
  }
  return(matrix(c("p_value","Conclusión",pval,h),2,2))
}
Prueba_homocedasticidad(RunSize,RunTime,Y_reg)

##      [,1]         [,2]               
## [1,] "p_value"    "0.428518409171384"
## [2,] "Conclusión" "No se rechaza Ho"

Hipótesis de Normalidad

\[H_0: \text{El modelo sigue una distribución Normal}\] \[H_a: \text{El modelo no sigue una distribución Normal}\]

Prueba_normalidad<-function(X,Y,regresion_lineal)
{
  residuos <- function(X,Y,regresion_lineal){
    yreg<- regresion_lineal(X)
    return(Y-yreg)
  }
  Residuos<- residuos(RunSize,RunTime,Y_reg)
  p<-shapiro.test(Residuos) 
  library(nortest)
  b<-lillie.test(Residuos)
  pval<-p$p.value
  pval2<-p$p.value
  if (pval > 0.05){
    resultado<-TRUE
    h<-c("No se rechaza Ho.")
  }
  else{
    resultado<-FALSE
    h<-c("Se rechaza Ho.")
  }
  return(matrix(c("p_value Shapiro",pval,"p_value Lillie",pval2,"Conclusión",h),2,3))
}
Prueba_normalidad(RunSize,RunTime,Y_reg)

##      [,1]                [,2]                [,3]               
## [1,] "p_value Shapiro"   "p_value Lillie"    "Conclusión"       
## [2,] "0.891701087186422" "0.891701087186422" "No se rechaza Ho."

Hipótesis de Independencia

\[H_0: \rho=0\] \[H_a: \rho\not= 0\]