### definir directorio
#setwd("~/AESpecialist/Regresion")
mtcars
##                      mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
## Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
## Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
## Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
## Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
## Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
## Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
## Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
## Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
## Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
## Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
## AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
## Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
## Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
## Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
## Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
## Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
## Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
## Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
## Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2
library(readxl)
library(MASS)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(nortest)
library(MVN)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
## sROC 0.1-2 loaded
attach(mtcars)## exploracion y resumen datos 
head(mtcars)
##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1
summary(mtcars)
##       mpg             cyl             disp             hp       
##  Min.   :10.40   Min.   :4.000   Min.   : 71.1   Min.   : 52.0  
##  1st Qu.:15.43   1st Qu.:4.000   1st Qu.:120.8   1st Qu.: 96.5  
##  Median :19.20   Median :6.000   Median :196.3   Median :123.0  
##  Mean   :20.09   Mean   :6.188   Mean   :230.7   Mean   :146.7  
##  3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.:326.0   3rd Qu.:180.0  
##  Max.   :33.90   Max.   :8.000   Max.   :472.0   Max.   :335.0  
##       drat             wt             qsec             vs        
##  Min.   :2.760   Min.   :1.513   Min.   :14.50   Min.   :0.0000  
##  1st Qu.:3.080   1st Qu.:2.581   1st Qu.:16.89   1st Qu.:0.0000  
##  Median :3.695   Median :3.325   Median :17.71   Median :0.0000  
##  Mean   :3.597   Mean   :3.217   Mean   :17.85   Mean   :0.4375  
##  3rd Qu.:3.920   3rd Qu.:3.610   3rd Qu.:18.90   3rd Qu.:1.0000  
##  Max.   :4.930   Max.   :5.424   Max.   :22.90   Max.   :1.0000  
##        am              gear            carb      
##  Min.   :0.0000   Min.   :3.000   Min.   :1.000  
##  1st Qu.:0.0000   1st Qu.:3.000   1st Qu.:2.000  
##  Median :0.0000   Median :4.000   Median :2.000  
##  Mean   :0.4062   Mean   :3.688   Mean   :2.812  
##  3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:4.000  
##  Max.   :1.0000   Max.   :5.000   Max.   :8.000
dim(mtcars)
## [1] 32 11
  1. Realizaremos test de normalidad bivariada

realizamos prueba de Mardia \[ H_0: \text{elvector aleatorio distribuye normal bivariado }\\H_a: \text{el vector aleatorio No distribuye normal bivariado}\]

mvn(mtcars)
## $multivariateNormality
##              Test          Statistic              p value Result
## 1 Mardia Skewness   394.609460003421 2.11444606575459e-05     NO
## 2 Mardia Kurtosis 0.0383039985304523    0.969445302761312    YES
## 3             MVN               <NA>                 <NA>     NO
## 
## $univariateNormality
##            Test  Variable Statistic   p value Normality
## 1  Shapiro-Wilk    mpg       0.9476  0.1229      YES   
## 2  Shapiro-Wilk    cyl       0.7533  <0.001      NO    
## 3  Shapiro-Wilk   disp       0.9200  0.0208      NO    
## 4  Shapiro-Wilk    hp        0.9334  0.0488      NO    
## 5  Shapiro-Wilk   drat       0.9459  0.1101      YES   
## 6  Shapiro-Wilk    wt        0.9433  0.0927      YES   
## 7  Shapiro-Wilk   qsec       0.9733  0.5935      YES   
## 8  Shapiro-Wilk    vs        0.6323  <0.001      NO    
## 9  Shapiro-Wilk    am        0.6251  <0.001      NO    
## 10 Shapiro-Wilk   gear       0.7728  <0.001      NO    
## 11 Shapiro-Wilk   carb       0.8511   4e-04      NO    
## 
## $Descriptives
##       n       Mean     Std.Dev  Median    Min     Max      25th   75th
## mpg  32  20.090625   6.0269481  19.200 10.400  33.900  15.42500  22.80
## cyl  32   6.187500   1.7859216   6.000  4.000   8.000   4.00000   8.00
## disp 32 230.721875 123.9386938 196.300 71.100 472.000 120.82500 326.00
## hp   32 146.687500  68.5628685 123.000 52.000 335.000  96.50000 180.00
## drat 32   3.596563   0.5346787   3.695  2.760   4.930   3.08000   3.92
## wt   32   3.217250   0.9784574   3.325  1.513   5.424   2.58125   3.61
## qsec 32  17.848750   1.7869432  17.710 14.500  22.900  16.89250  18.90
## vs   32   0.437500   0.5040161   0.000  0.000   1.000   0.00000   1.00
## am   32   0.406250   0.4989909   0.000  0.000   1.000   0.00000   1.00
## gear 32   3.687500   0.7378041   4.000  3.000   5.000   3.00000   4.00
## carb 32   2.812500   1.6152000   2.000  1.000   8.000   2.00000   4.00
##            Skew    Kurtosis
## mpg   0.6106550 -0.37276603
## cyl  -0.1746119 -1.76211977
## disp  0.3816570 -1.20721195
## hp    0.7260237 -0.13555112
## drat  0.2659039 -0.71470062
## wt    0.4231465 -0.02271075
## qsec  0.3690453  0.33511422
## vs    0.2402577 -2.00193762
## am    0.3640159 -1.92474143
## gear  0.5288545 -1.06975068
## carb  1.0508738  1.25704307
cor.test(disp,wt, data="mtcars")## coeficiente pearson 
## 
##  Pearson's product-moment correlation
## 
## data:  disp and wt
## t = 10.576, df = 30, p-value = 1.222e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7811586 0.9442902
## sample estimates:
##       cor 
## 0.8879799

Para este test la hipotesis sera la siguiente \[ H_0: \rho= 0\\H_a:\rho \not= 0\] EL test de mardia confirma la normalidad de los datos de mt cars pero esta no es normalidad bivariada, el test confirma normalidad univariada ya que el skewness no se cumple en los datos mtcars. La correlacion de pearson nos da valor de 0.8, la corrimos sin embargo ya que la prueba de mardia resulto en univariada el coeficiente de correlacion mas apropiado es el de spearman, no es necesario interpretar.

La decisión se toma a partir del valor p, en este caso con un alfa de 0.01

  1. Realizamos test de normalidad bivariada de spearman \[ H_0: \rho= 0\\H_a:\rho \not= 0\] dado que el test de mardia del punto 1 dio que el vector no es normal bivariado, el coeficiente de spearman no parametrico es el mas apropiado para evaluar la correlacion de las variables
cor.test(disp, wt, data="mtcars", method = "spearman")## coeficiente spearman 
## Warning in cor.test.default(disp, wt, data = "mtcars", method = "spearman"):
## Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  disp and wt
## S = 558.11, p-value = 3.346e-12
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8977064

La decisión se toma a partir del valor p, en este caso con un alfa de 0.01, se rechaza H_0 El coeficiente de spearman da 0.9, puedo decir qu existe una relacion positiva y significativa entre dist y wt a medida qye aumenta una aumenta la otra.

  1. Para este punto adjuntaremos nuestr base de babies
babies<-read_excel("C:\\Users\\yarvis\\Documents\\AESpecialist\\Regresion\\babys.xlsx")#llamar datos 
attach(babies)
head(babies)
## # A tibble: 6 x 5
##   Talla  Edad Talla_nacer Peso_nacer Torax_nacer
##   <dbl> <dbl>       <dbl>      <dbl>       <dbl>
## 1  57.5    78        48.2       2.75        29.5
## 2  52.8    69        45.5       2.15        26.3
## 3  61.3    77        46.3       4.41        32.2
## 4  67      88        49         5.52        36.5
## 5  53.5    67        43         3.21        27.2
## 6  62.7    80        48         4.32        27.7
summary(babies)
##      Talla            Edad      Talla_nacer      Peso_nacer     Torax_nacer   
##  Min.   :52.80   Min.   : 67   Min.   :43.00   Min.   :2.150   Min.   :26.30  
##  1st Qu.:56.20   1st Qu.: 74   1st Qu.:46.30   1st Qu.:2.750   1st Qu.:27.70  
##  Median :61.30   Median : 78   Median :48.00   Median :3.710   Median :28.70  
##  Mean   :60.97   Mean   : 81   Mean   :48.78   Mean   :3.631   Mean   :29.63  
##  3rd Qu.:67.00   3rd Qu.: 88   3rd Qu.:49.00   3rd Qu.:4.320   3rd Qu.:30.30  
##  Max.   :69.20   Max.   :102   Max.   :58.00   Max.   :5.520   Max.   :36.50
dim(babies)
## [1] 9 5
summary(babies)
##      Talla            Edad      Talla_nacer      Peso_nacer     Torax_nacer   
##  Min.   :52.80   Min.   : 67   Min.   :43.00   Min.   :2.150   Min.   :26.30  
##  1st Qu.:56.20   1st Qu.: 74   1st Qu.:46.30   1st Qu.:2.750   1st Qu.:27.70  
##  Median :61.30   Median : 78   Median :48.00   Median :3.710   Median :28.70  
##  Mean   :60.97   Mean   : 81   Mean   :48.78   Mean   :3.631   Mean   :29.63  
##  3rd Qu.:67.00   3rd Qu.: 88   3rd Qu.:49.00   3rd Qu.:4.320   3rd Qu.:30.30  
##  Max.   :69.20   Max.   :102   Max.   :58.00   Max.   :5.520   Max.   :36.50
boxplot(babies)

segun la estrategia revisada en clase realizaremos los siguientes pasos 1. Definimos modelo teórico \[ Y_i= \beta_0+\beta_1X_i+\varepsilon_i\] probamos si existe un modelo RLS

\[ H_0: \beta_1= 0\\H_a:\beta_1 \not= 0\] Si beta es cero tendremos un modleo nulo y no habria varable regresora

baby.fit11<-lm(Talla_nacer~Torax_nacer)
baby.fit11
## 
## Call:
## lm(formula = Talla_nacer ~ Torax_nacer)
## 
## Coefficients:
## (Intercept)  Torax_nacer  
##     42.2965       0.2187

Estamos ante un modelo RLS 2. Definir supuestos para los errores

residuals(baby.fit11)
##          1          2          3          4          5          6          7 
## -0.5486158 -2.5487274 -3.0391466 -1.2796215 -5.2455710 -0.3549286 -0.4861576 
##          8          9 
##  4.0764122  9.4263563
qf(0.99,1,7)# prueba f 
## [1] 12.24638

en la prueba f tomamos los grados de libertad de n-2 para el denominador en el casode nuestros datos seria 7

anova(baby.fit11)
## Analysis of Variance Table
## 
## Response: Talla_nacer
##             Df  Sum Sq Mean Sq F value Pr(>F)
## Torax_nacer  1   3.713  3.7131  0.1721 0.6907
## Residuals    7 151.023 21.5746

EL valor de f critico es mayor que el valor de 12.24 por tanto podemos decir que la variable torax es regresora para talla nacer

  1. Estimar los coeficientes del módelo \[\beta_0 , \beta_1\]
baby.fit11  
## 
## Call:
## lm(formula = Talla_nacer ~ Torax_nacer)
## 
## Coefficients:
## (Intercept)  Torax_nacer  
##     42.2965       0.2187

Solo estimando los coeficientes del modelo podemos saber como son los residuales

  1. Expresar el modelo ajustado

\[ \hat y= \hat \beta_0+\hat\beta_1X_i\]

  1. Generar los n residuales del modelo \[\hat\varepsilon_i= y_i-\hat y_i\]
resid.torax<-residuals(baby.fit11)  ### Residuales (errores estimados)
resid.torax
##          1          2          3          4          5          6          7 
## -0.5486158 -2.5487274 -3.0391466 -1.2796215 -5.2455710 -0.3549286 -0.4861576 
##          8          9 
##  4.0764122  9.4263563

La variabilidad del modelo no dice si tenemos un modelo de regresion simple L suma de cuadrados del error es mayor que la suma de cuadrados del modelo Estimación de los errores

  1. Con los residuales ,validar los supuestos sobre los errores

¿Provienen de una distibucion normal?

\[H_0:\text{los residuales provienen de una distribucion normal}\\H_a:\text{los residuales no provienen de una diistribucion normal}\]

shapiro.test(residuals(baby.fit11)) ## normalidad 
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(baby.fit11)
## W = 0.86878, p-value = 0.1193

Significancia de la prueba 0.01

mediante el valor de p decidimos que rechazamos y no encontramos normalidad en los datos

lillie.test(residuals(baby.fit11)) ### normalidad
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuals(baby.fit11)
## D = 0.31033, p-value = 0.01258

Significancia de la prueba 0.01

si no fuese normal no se puede construir prueba f

mediante el valor de p decidimos que rechazamos y no encontramos normalidad en los datos

Varianza constante de los errores

\[H_0:\text{los residuales tienen varianza constante}\\H_a:\text{los residuales no tienen varianza constante}\]

bptest(Talla_nacer~ Torax_nacer) ### Homocedasticidad Breusch-Pagan 
## 
##  studentized Breusch-Pagan test
## 
## data:  Talla_nacer ~ Torax_nacer
## BP = 0.28262, df = 1, p-value = 0.595

Significancia de la prueba 0.01 Se rechaza H_0 podemos concluir que los residuales no tiene varianza constante

gqtest(lm(Talla_nacer~ Torax_nacer))
## 
##  Goldfeld-Quandt test
## 
## data:  lm(Talla_nacer ~ Torax_nacer)
## GQ = 12.751, df1 = 3, df2 = 2, p-value = 0.0736
## alternative hypothesis: variance increases from segment 1 to 2

Significancia de la prueba 0.01

Se rechaza H_0 podemos concluir que los residuales no tiene varianza constante

Independencia de los errores La hipotesis para independencia se toma de la siguiente manera,

\[H_0:\text{los errores son independientes}\\H_a:\text{los errores no son independientes}\] aun que en este test de independencia la Ha deberia ser que la autocorrelacion es mayor que cero, por eso es mas factible asumir las anteriores hipotesis. Ya que si tengo variables independientes podemos asumir correlación cero.

dwtest(Talla_nacer ~ residuals(baby.fit11)) # Durbin Watson independencia
## 
##  Durbin-Watson test
## 
## data:  Talla_nacer ~ residuals(baby.fit11)
## DW = 2.0253, p-value = 0.3787
## alternative hypothesis: true autocorrelation is greater than 0

para construir distribucion chi_cuadrado es por que las distribuciones de los errores son independientes

Significancia de l prueba 0.05

Si escojo el peor modelo me comprometo a seguir buscando un modelo

  1. Si se cumplen los supuestos interpretar los parametros del modelo ajustado \[Y= \beta_0+\beta_1*Toraxnacer\]
anova(baby.fit11)
## Analysis of Variance Table
## 
## Response: Talla_nacer
##             Df  Sum Sq Mean Sq F value Pr(>F)
## Torax_nacer  1   3.713  3.7131  0.1721 0.6907
## Residuals    7 151.023 21.5746
babies
## # A tibble: 9 x 5
##   Talla  Edad Talla_nacer Peso_nacer Torax_nacer
##   <dbl> <dbl>       <dbl>      <dbl>       <dbl>
## 1  57.5    78        48.2       2.75        29.5
## 2  52.8    69        45.5       2.15        26.3
## 3  61.3    77        46.3       4.41        32.2
## 4  67      88        49         5.52        36.5
## 5  53.5    67        43         3.21        27.2
## 6  62.7    80        48         4.32        27.7
## 7  56.2    74        48         2.31        28.3
## 8  68.5    94        53         4.3         30.3
## 9  69.2   102        58         3.71        28.7