Supuestos

1.-La relación entre Y y X es lineal.
2.-Las Xi son variables no estocásticas (no aleatorias) cuyos valores son fijos (o controlados por el experimentador).
3.-La media del término error es cero.
4.-La varianza del término error es igual en todos los niveles de X. Es el principio de homocedasticidad.
-5.-No existe autocorrelación entre el término error de diferentes niveles de X.(con durbin watason
**6.-La covarianza entre el término error y la variable explicativa X es cero.
7.-El número de observaciones (el tamaño de la muestra) debe ser superior al número de parámetros a estimar.
8.-Los valores de X en una muestra dada no pueden ser todos iguales, teóricamente la var(X) debe ser un número finito positivo.
*9.-El modelo está correctamente especificado.
-10.-No existe multicolinealidad perfecta. (para regresiones multiples )

Primer ejercicio

Importar Dataset

#llamar la biblioteca de datasets
library(help="datasets")
data("iris")
head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
#solo voy a tomar del 51 al 150 para eliminar setosas y se vea mejor la relación lineal
iris<-iris[51:150,]

Explorar los datos

#sacar la correlación entre las variables
cor(iris[,1:4])
##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length    1.0000000   0.5538548    0.8284787   0.5937094
## Sepal.Width     0.5538548   1.0000000    0.5198023   0.5662025
## Petal.Length    0.8284787   0.5198023    1.0000000   0.8233476
## Petal.Width     0.5937094   0.5662025    0.8233476   1.0000000

7.-El número de observaciones (el tamaño de la muestra) debe ser superior al número de parámetros a estimar.

se cumple por default, casi siempre se va a cumplir, n> número de betas, si quiero evaluar si variables debe haber al menos 6 observaciones, en este caso hay 100 observaciones para un factor 

str(iris)
## 'data.frame':    100 obs. of  5 variables:
##  $ Sepal.Length: num  7 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 ...
##  $ Sepal.Width : num  3.2 3.2 3.1 2.3 2.8 2.8 3.3 2.4 2.9 2.7 ...
##  $ Petal.Length: num  4.7 4.5 4.9 4 4.6 4.5 4.7 3.3 4.6 3.9 ...
##  $ Petal.Width : num  1.4 1.5 1.5 1.3 1.5 1.3 1.6 1 1.3 1.4 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 2 2 2 2 2 2 2 2 2 2 ...
1.-La relación entre Y y X es lineal.

Hay una clara relación lineal

8.-Los valores de X en una muestra dada no pueden ser todos iguales, teóricamente la var(X) debe ser un número finito positivo.

9.-El modelo está correctamente especificado.

teóricamente sí existe una relación entre el pétalo y el sépalo de una flor 

plot(y=iris$Sepal.Length, x=iris$Petal.Length ,
     col="lightblue", 
     pch=18,
     main="Relación Sépalo Pétalo",
     xlab="Longitud del pétalo cm",
     ylab="Longitud del sépalo cm",)


Manualmente alpha y beta

#desviación estandar en x
sd(iris$Petal.Length)
## [1] 0.8255785
#desviación estandar en y
sd(iris$Sepal.Length)
## [1] 0.6628344
correlacionXY<-cor(iris[,1:4])
print(correlacionXY)
##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length    1.0000000   0.5538548    0.8284787   0.5937094
## Sepal.Width     0.5538548   1.0000000    0.5198023   0.5662025
## Petal.Length    0.8284787   0.5198023    1.0000000   0.8233476
## Petal.Width     0.5937094   0.5662025    0.8233476   1.0000000
correlacionXY<-cor(iris[,1:4])[3,1] #la posición columna 3,fila 1
print(correlacionXY)
## [1] 0.8284787
#valor de la pendiente
correlacionXY<-cor(iris[,1:4])[3,1]
desvX<-sd(iris$Petal.Length)
desvY<-sd(iris$Sepal.Length)

pendiente<-correlacionXY*(desvY/desvX)
print(pendiente)
## [1] 0.6651629
#calcular el valor de la pendiente 2
covarianzaXY<-cov(iris$Petal.Length,iris$Sepal.Length)
varianzaX<-var(iris$Petal.Length)

pendiente2<-covarianzaXY/varianzaX
print(pendiente2)
## [1] 0.6651629
#cálculo del intercepto
mediaY<-mean(iris$Sepal.Length)
mediaX<-mean(iris$Petal.Length)

intercepto<-mediaY-(pendiente2*mediaX)
print(intercepto)
## [1] 2.998711

Hacer el modelo

modelo<-lm(Sepal.Length~Petal.Length,iris)
#ver los interceptos
modelo
## 
## Call:
## lm(formula = Sepal.Length ~ Petal.Length, data = iris)
## 
## Coefficients:
##  (Intercept)  Petal.Length  
##       2.9987        0.6652

Ver que el p valor es menor que 0.05 lo que significa que si hay una pendiente

summary(modelo)
## 
## Call:
## lm(formula = Sepal.Length ~ Petal.Length, data = iris)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.09194 -0.26570  0.00761  0.21902  0.87502 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.99871    0.22593   13.27   <2e-16 ***
## Petal.Length  0.66516    0.04542   14.64   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3731 on 98 degrees of freedom
## Multiple R-squared:  0.6864, Adjusted R-squared:  0.6832 
## F-statistic: 214.5 on 1 and 98 DF,  p-value: < 2.2e-16

Análizar los residuos

residuos=residuals(modelo)
Media del error debe ser cero

3.-La media del término error es cero. 

mean(residuos)
## [1] -3.057179e-17

Los errores son independientes. los errores se distribuyen en forma normal

par(mfrow=c(2,1))
hist(residuos, col="yellow")
boxplot(residuos, bty="l", range=1.5, col="yellow", horizontal=T,xlab="residuos")

par(mfrow=c(1,1))
hist(residuos,
     breaks=25,
     col="lightblue",
     )

hist(residuos,
     breaks=25,
     col="lightblue",
     freq=FALSE,
     )

lines(density(residuos), lwd = 2, col = 'blue')

hist(residuos,
     breaks=25,
     col="lightblue",
     freq=FALSE,
     )

lines(density(residuos), lwd = 2, col = 'blue')

#tendencia de curva normal
x <- seq(min(residuos), max(residuos), length = 40)
f <- dnorm(x, mean = mean(residuos), sd = sd(residuos))
lines(x, f, col = "red", lwd = 2)


Test de normalidad h0= distribución normal, se busca valores mayores a 0.05

shapiro.test(residuos) #menos de 50 observaciones
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.99268, p-value = 0.8682
#para más de 50 observaciones
#kolmorogov
ks.test(residuos,"pnorm", mean(residuos), sd(residuos))
## Warning in ks.test(residuos, "pnorm", mean(residuos), sd(residuos)): ties should
## not be present for the Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.043847, p-value = 0.9907
## alternative hypothesis: two-sided
ks.test(residuos, "pnorm") #creo que es mejor la de arriba
## Warning in ks.test(residuos, "pnorm"): ties should not be present for the
## Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  residuos
## D = 0.23614, p-value = 2.868e-05
## alternative hypothesis: two-sided
# Alternativamente conviene usar la modificación de Lilliefors a este test.
# Esta corrección considera que los parámetros son estimados, a diferencia
# del «ks» «a secas»:

library(nortest)
lillie.test(residuos)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuos
## D = 0.043847, p-value = 0.9076
# realizar la prueba Anderson-Darling para comprobar la normalidad
ad.test (residuos)
## 
##  Anderson-Darling normality test
## 
## data:  residuos
## A = 0.19852, p-value = 0.8841

Prueba de homocedasticidad

4.-La varianza del término error es igual en todos los niveles de X. Es el principio de homocedasticidad.

Se busca un pvalor mayor a 0.05 ya que la h0=las pruebas son homogeneas

#install.packages('lmtest')
library(lmtest)
bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 0.04988, df = 1, p-value = 0.8233
predichos=fitted(modelo)

plot(residuos~predichos)

plot(residuos)

#grafico unido
par(mfcol = c(2, 2))
plot(modelo)

par(mfcol = c(1, 1))
plot(modelo,1)

6.-La covarianza entre el término error y la variable explicativa X es cero.

cov(residuos,iris$Petal.Length)
## [1] -3.589982e-17
cor(residuos,iris$Petal.Length)
## [1] -1.171453e-16
plot(residuos,iris$Petal.Length)