1. Preparacion

En esta sección se muestra las librerías a utilizar y funciones que se requieren para realizar las simulaciones

a. Lista de paquetes a utilizar

library(Matrix)
library(stats)
library(broom)
library(dplyr)
library(ggfortify)
library(MASS)
library(EnvStats)
library(lmtest)
library(DescTools)
library(gvlma)
library(psych)
library(knitr)
library(Metrics)

b. Función Generadora de Datos: cumple los supuestos

Se procede a generar una función R que cumple los supuestos del modelo de regresión lineal:

\[E^{\theta}(y_{i}|z_{i}) = EL^{\theta}(y_{i}|z_{i})= \beta^{\prime}z_{i}= \displaystyle\sum_{j=1}^{q}\beta _{j}z_{ij} \]

# Proceso Generador de datos
f_ProcesoGeneradorDatos = function(pN,pB0,pB1,pB2,pB3,pSigma){
  psd = rep(pSigma,pN)
    
  Z0 = rep(1,pN)
  
  Z1 = runif(pN,20,80)
  
  Z2 = rexp(pN,0.01)
  
  Z3 = rnorm(pN,mean=40,sd=5)
  
  Y = rnorm(pN,mean = (pB0 + pB1*Z1 + pB2*Z2 + pB3*Z3),sd= psd)
  
  M = as.matrix(cbind (Y,Z0,Z1,Z2,Z3))
  
  return(M)
}

Esto nos permite simular datos como por ejemplo: E(Y/Z)~ N(10Z0 + 0.3Z1 + 0.4Z2 - 1.2Z3, 8)

Parámetros :

B0     = 10
B1     = 0.3
B2     = 0.4
B3     = 1.2
Sigma  = 8

Mostramos las primeras 20 observaciones generadas por la función:

M = f_ProcesoGeneradorDatos(20,B0,B1,B2,B3,Sigma)
M
##               Y Z0       Z1         Z2       Z3
##  [1,]  58.61575  1 53.37349   9.048203 36.34459
##  [2,] 123.50782  1 20.99971 134.921563 44.24281
##  [3,] 162.10263  1 38.21572 219.637337 37.65139
##  [4,]  99.36736  1 36.88638 102.085008 43.74056
##  [5,] 104.00104  1 41.81047  25.357506 46.17104
##  [6,]  98.97224  1 48.01159  29.248950 53.74811
##  [7,] 107.33137  1 78.06965  80.426349 35.07459
##  [8,]  74.72967  1 31.46951  23.631980 41.83440
##  [9,] 131.78014  1 27.11963 150.972006 40.69657
## [10,] 164.51520  1 71.30100 157.845088 45.78310
## [11,]  80.46052  1 75.75088   1.113151 34.89797
## [12,] 112.30219  1 50.24977  99.328566 37.24104
## [13,]  98.66661  1 53.29017  20.422773 46.02805
## [14,]  93.36489  1 50.99708  71.492748 32.47757
## [15,] 150.16766  1 58.89884 153.107629 41.30495
## [16,] 132.15959  1 26.60537 188.338647 39.56592
## [17,] 209.73811  1 42.67583 294.258413 49.38902
## [18,]  56.14117  1 39.03198  22.531373 37.36127
## [19,]  96.85721  1 25.09469 108.079396 37.13552
## [20,]  71.87900  1 27.33701   7.499001 36.98669

Mostramos un gráfico para observar la distribución de las variables

multi.hist(x = M[,-2])

c. Función Generadora de Datos: no cumple los supuestos

También creamos una segunda función que no cumple los supuestos de un modelo de regresión lineal:

# Proceso Generador de datos: No cumple
f_ProcesoGeneradorDatos2 = function(pN,pB0,pB1,pB2,pB3,pSigma){
  psd = seq(0, 10, length.out = pN)  + pSigma
  
  Z0 = rep(1,pN)
  
  Z1 = runif(pN,20,80)
  
  # Correlación
  Z2 = rexp(pN,0.01) + Z1
  
  Z3 = rnorm(pN,mean=40,sd=5)
  
  # No normalidad
  Y = pB0 + sqrt(pB1)*Z1 + (pB2^2)*Z2 + pB3*Z3 + rexp(pN, psd/1000)
  
  M = as.matrix(cbind (Y,Z0,Z1,Z2,Z3))
  
  return(M)
}

Esta función genera datos que no cumplen con los supuestos de un modelo de regresión lineal

M = f_ProcesoGeneradorDatos2(20,B0,B1,B2,B3,Sigma)
M
##               Y Z0       Z1        Z2       Z3
##  [1,] 245.43147  1 46.47229 222.14955 37.08346
##  [2,] 170.49893  1 42.62876 491.46441 41.25707
##  [3,] 143.82893  1 28.92184  34.46972 40.37954
##  [4,] 142.99672  1 22.41173  30.31685 40.20813
##  [5,] 305.62195  1 43.31085  60.64287 29.97321
##  [6,] 160.80544  1 66.61025 390.24539 40.99429
##  [7,] 153.54871  1 74.80918  93.10748 50.02064
##  [8,] 369.36865  1 31.23908 207.73259 41.62709
##  [9,] 115.90900  1 44.49722  60.29334 43.03583
## [10,] 175.02683  1 20.81276 230.10312 41.61067
## [11,] 158.39839  1 28.50718  30.15324 45.33033
## [12,] 208.75399  1 64.87750  75.91679 44.75717
## [13,] 134.96707  1 20.03465 170.94656 41.89297
## [14,] 104.02692  1 47.79180  98.60955 37.31221
## [15,]  92.75251  1 50.88560  89.89255 33.18095
## [16,] 167.65504  1 50.66531 187.76888 29.49964
## [17,] 158.10667  1 44.11114 346.02237 46.75996
## [18,] 137.49537  1 59.48950  87.94536 38.18702
## [19,] 129.62073  1 64.71303 178.18444 42.01113
## [20,] 174.15530  1 54.32674 199.42512 47.26528
multi.hist(x = M[,-2])

d. Función que estima los parámetros de modelo lineal

Esta función estima los parámetros: \(\beta_{0},\phantom{a}\beta_{1} ,\phantom{a}\beta_{2},\phantom{a} \beta_{3}\phantom{a} y \phantom{a}\sigma\)

f_RegresionLineal= function(pM){
  
  Y  = pM[,1]
  Z  = pM[,2:5]
  N  = nrow(Z)
  k  = ncol(Z)
  
  B_hat = solve(t(Z)%*%Z)%*%(t(Z)%*%Y)
  
  Y_hat = Z%*%B_hat
  
  e_hat = Y - Y_hat
  
  SCR = t(e_hat)%*%e_hat
  sCT = t(Y)%*%Y- N*(mean(Y)^2)
  
  sigma_mco = sqrt(SCR/(N - k))
  
  sigma_mv  = sqrt(SCR/(N))
  
  xBO    = B_hat[1,1]
  xB1    = B_hat[2,1]
  xB2    = B_hat[3,1]
  xB3    = B_hat[4,1]
  xs_mco = sigma_mco
  xs_mv  = sigma_mv
  XR2    = 1 - SCR/sCT
  
  
  B =  data.frame( BO = xBO , B1 = xB1  ,B2 = xB2   ,B3 = xB3 , S_MCO = xs_mco, S_MV = xs_mv ,R2 = XR2,row.names = NULL )
  
  return(B)
  
}

Con ellos podemos estimar los parámetros del modelo dado los parámetros poblacionales

B0     = 10
B1     = 0.3
B2     = 0.4
B3     = 1.2
Sigma  = 8


M = f_ProcesoGeneradorDatos(200,B0,B1,B2,B3,Sigma)

P = f_RegresionLineal(M)
P

2. Supuestos de un modelo de regresión lineal

Para simular la validación de los supuestos, generamos 100 observaciones que cumplen con 4 supuestos de un modelo lineal

B0     = 10
B1     = 0.3
B2     = 0.4
B3     = 1.2
Sigma  = 8


M = f_ProcesoGeneradorDatos(200,B0,B1,B2,B3,Sigma)

mrl  = lm(Y~Z1+Z2+Z3,data = as.data.frame(M))

coefficients(mrl)
## (Intercept)          Z1          Z2          Z3 
##  11.0620101   0.3701391   0.4055110   1.0771049

Para validar los supuestos es importante analizar los residuos que se generan, es decir, la diferencia entre la valor real de la variable respuesta y el valor estimado de la misma:

\[u_{i} = \widehat{y_{i}} - y_{i}\]

residuos = residuals(mrl)


head(residuals(mrl), 10)
##          1          2          3          4          5          6          7 
## -0.9961438 -2.4953599  3.8917899  8.3467562 -5.5513838 -0.9301819 -2.1572065 
##          8          9         10 
## -0.1383001 -2.7033042  0.5462970

a. Supuesto A1

El primers supuesto es el que corresponde a la linealidad del modelo:

\[E^{\theta}(y_{i}|z_{i}) = EL^{\theta}(y_{i}|z_{i})= \beta^{\prime}z_{i}= \displaystyle\sum_{j=1}^{q}\beta _{j}z_{ij} \]

Se deduce que este supuesto puede ser verificado si la media de los residuos es 0

\[E^{\theta}(u_{i}|z_{i}) = 0\]

mean(residuos)
## [1] -2.489849e-16
plot(residuos) + abline (h = 0, col = 'red', lwd = 2)

## integer(0)
autoplot(mrl,1)

Se cumple el supuesto A1

b. Supuesto A2

El segundo supuesto indica que la matriz de datos \(Z\) debe tener rango completo:

\[Rango(Z) = q \;\;(n > q)\]

O equivalentemente :

\[Rango(Z'Z)=q,; det(Z'Z)\neq 0 ; Z'Z\phantom{a}invertible\]

Z = M[,2:5]
rankMatrix(Z)[1]
## [1] 4
det(t(Z)%*%Z)
## [1] 7.530774e+16
solve(t(Z)%*%Z)
##               Z0            Z1            Z2            Z3
## Z0  4.420461e-01 -9.254413e-04 -8.151921e-05 -9.569006e-03
## Z1 -9.254413e-04  1.799860e-05  1.313269e-07  7.670032e-07
## Z2 -8.151921e-05  1.313269e-07  6.264550e-07  4.665214e-07
## Z3 -9.569006e-03  7.670032e-07  4.665214e-07  2.362716e-04

Se cumple el supuesto A2

c. Supuesto A3

El tercer supuesto indica la homocedasticidad de la variable respuesta \(y_{i}\) condicionada a los valores de \(z_{i}\) debe ser constante

\[Var^{\theta}(y_{i}|z_{i})=\sigma^{2}\phantom{abc}para\phantom{a}todo\phantom{abc}i=1,\ldots,n.\]

Este supuesto se puede verificar si se cumplen estas dos condiciones:

  1. Homocedasticidad en los residuos : \(E^{\theta}(u_{i}u_{i}|z_{i}) = {\sigma}^2\)
  2. No existe correlación en los residuos : \(E^{\theta}(u_{i}u_{j}|z_{i}) = 0\)
# A3.1 Homocedasticidad de los residuos Var(ui) = cte
# studentized Breusch-Pagan test

#Ho : los residuos tienen varianza constante.
#H1 : los residuos no tienen varianza constante.
#p-valor es mayor de 0.05, no se rechaza HO
bptest(mrl)
## 
##  studentized Breusch-Pagan test
## 
## data:  mrl
## BP = 0.61881, df = 3, p-value = 0.8921
autoplot(mrl,3)

Se verifica que los residuos son homocedásticos

# A3.1 Autorrelacion entre los residuos : Cov(ui,uj) = 0

#Ho : Los residuos no estan autocorrelacionados
#Ha:  Los residuos están autocorrelacionados
#p-valor es mayor de 0.05, no se rechaza HO

DurbinWatsonTest(mrl)
## 
##  Durbin-Watson test
## 
## data:  mrl
## DW = 2.0568, p-value = 0.6546
## alternative hypothesis: true autocorrelation is greater than 0
# Todos los barras se encuetra entre de las bandas
acf(residuos )

Se verifica que los residuos no se encuentran correlacionados

d. Supuesto A4

Este supuesto corresponde a normalidad condicionada de la variable de respuesta

\[ y_{i}|z_{i} \sim i.i.N(\beta^{\prime}z_{i},\sigma^{2}) \]

Este supuesto es equivalente a la normalidad de los residuos:

\[ u_{i}|z_{i} \sim i.i.N(0,\sigma^{2}) \]

#Kolmogorov-Smirnov test

#H0: los residuos proceden de una distribución normal
#H1: los residuos no proceden de una distribución normal
#p-valor es mayor de 0.05, no se rechaza HO

ks.test(residuals(mrl),pnorm (0,Sigma))
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  residuals(mrl) and pnorm(0, Sigma)
## D = 0.52, p-value = 0.9652
## alternative hypothesis: two-sided
# Graficos
hist( residuos)

qqPlot(residuos) 

autoplot(mrl,2)

Se verfica que los residuos tienen uan distribucion normal como media 0 y varianza constante

Podemos resumir la verificación de los supuetsos A1, A2, A3 y A4 en el siguiente gráfico:

# Graficos para verificar los principales supuestos
#plot(gvlma(mrl))

3. Propiedades de los Estimadores

En esta sección se verifican algunas propiedades de los estimadores

a. Consistencia

# P1 CONSISTENCIA plim(B2_hat) = B2
N = 10000


M = f_ProcesoGeneradorDatos(N,B0,B1,B2,B3,Sigma)

R = data.frame()


for(i in seq(10, N, by = 1)){

  P = cbind(f_RegresionLineal(M[1:i,]),"N"=i)

  R = rbind(R,P)
 
}


ggplot(data=R, aes(x=N, y=B2)) +
  geom_line(color="blue") + geom_hline(yintercept = B2, color = "red")

b. Insesgadez

# P2 INSESGADEZ E(B2_hat) = B2, aplicamos TCL
N = 50
n = 1000

R = data.frame()


for(i in seq(1, n, by = 1)){
  
  M = f_ProcesoGeneradorDatos(N,B0,B1,B2,B3,Sigma)
  
  P = cbind(f_RegresionLineal(M),"N"= N)
  
  R = rbind(R,P)

}

hist(R$B2,freq = FALSE,breaks =100, xlab = sprintf("B= %f y B_hat = %f",B2,mean(R$B2)))  
lines(density(R$B2), lwd = 2, col = 'red')
abline(v = mean(R$B2) ,lwd = 2, col = 'blue')

c. La varianza estimada por MV es sesgada y la de MCO es insesgada

N = 50
n = 1000

R = data.frame()


for(i in seq(1, n, by = 1)){
  
  M = f_ProcesoGeneradorDatos(N,B0,B1,B2,B3,Sigma)
  
  P = cbind(f_RegresionLineal(M),"N"= N)
  
  R = rbind(R,P)
  
}

hist(R$S_MCO,freq = FALSE,breaks =100, xlab = sprintf("Sigma= %f y S_MCO = %f",Sigma,mean(R$S_MCO)))  
abline(v = Sigma ,lwd = 2, col = 'blue')
abline(v = mean(R$S_MCO) ,lwd = 2, col = 'green')

hist(R$S_MV,freq = FALSE,breaks =100, xlab = sprintf("Sigma= %f y S_MV = %f",Sigma,mean(R$S_MV)))  
abline(v = Sigma ,lwd = 2, col = 'blue')
abline(v = mean(R$S_MV) ,lwd = 2, col = 'green')

4. Predicción de un modelo que cumple los supuestos

En esta sección se hace uso de las propiedades asintóticas de los estimadores

B0     = 10
B1     = 0.3
B2     = 0.4
B3     = 1.2
Sigma  = 8


# Datos de la muestra (entrenamiento)
M1 = f_ProcesoGeneradorDatos(70,B0,B1,B2,B3,Sigma)

# Datos fuera de la muestra (test)
M2 = f_ProcesoGeneradorDatos(30,B0,B1,B2,B3,Sigma)

# Regresion Lineal
mrl  = lm(Y~Z1+Z2+Z3,data = as.data.frame(M1))

summary(mrl)
## 
## Call:
## lm(formula = Y ~ Z1 + Z2 + Z3, data = as.data.frame(M1))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.314  -5.878   1.275   4.969  13.676 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.671455   7.259393   0.368    0.714    
## Z1          0.265144   0.051815   5.117 2.89e-06 ***
## Z2          0.395661   0.006466  61.192  < 2e-16 ***
## Z3          1.419173   0.186131   7.625 1.24e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.944 on 66 degrees of freedom
## Multiple R-squared:  0.9849, Adjusted R-squared:  0.9842 
## F-statistic:  1433 on 3 and 66 DF,  p-value: < 2.2e-16
confint(mrl)
##                   2.5 %     97.5 %
## (Intercept) -11.8223914 17.1653004
## Z1            0.1616919  0.3685966
## Z2            0.3827515  0.4085705
## Z3            1.0475515  1.7907953

Las predicción dentro y fuera de la muestra

# Prediccion dentro de la muestra
F = fitted(mrl,interval)
M1 = cbind(M1, Y_hat = F )

head(M1, n = 10)
##            Y Z0       Z1        Z2       Z3     Y_hat
## 1   86.63781  1 55.19162  15.60118 38.82631  78.57923
## 2  213.91993  1 45.56357 368.62192 34.45080 209.49335
## 3  101.20189  1 58.86051  46.29636 47.45368 103.94065
## 4  133.85926  1 51.70315 180.58724 39.33129 143.64950
## 5  110.81754  1 72.64153  85.74488 40.16805 112.86328
## 6   85.17545  1 25.73945  28.30841 41.56195  79.68028
## 7  109.82021  1 51.32938 140.09526 32.66965 118.07528
## 8   79.43354  1 56.53064  25.90246 38.87194  83.07484
## 9  103.26120  1 41.68332  83.75998 39.83708 103.39983
## 10  69.50813  1 75.41742  26.45760 33.30779  80.40572
# Prediccion fuera de la muestra
P = predict(mrl, as.data.frame(M2[,3:5]),interval = 'confidence')
M2 = cbind(M2, P )

head(M2, n = 10)
##            Y Z0       Z1         Z2       Z3       fit       lwr       upr
## 1  110.17360  1 54.62586  35.720166 39.55591  87.42497  85.40005  89.44988
## 2   74.71601  1 65.29735  21.567264 40.16439  85.51823  82.88170  88.15475
## 3   98.69709  1 69.90459  66.318432 43.92076 109.77705 106.87840 112.67571
## 4  140.00625  1 50.48341 145.274998 45.15306 137.61651 134.97526 140.25777
## 5  127.26910  1 72.09972  75.541622 49.26971 121.59942 117.45492 125.74391
## 6   77.44540  1 48.68116  41.404326 38.76597  86.97670  85.04608  88.90732
## 7   65.54845  1 35.19314   5.034263 47.16891  80.93544  76.84153  85.02936
## 8  181.20495  1 76.45148 259.013611 41.05208 183.68373 180.46440 186.90307
## 9  197.45009  1 56.92082 312.241382 43.29901 202.75422 199.59932 205.90913
## 10 119.56210  1 44.62702 134.711311 39.02293 123.18437 121.34925 125.01948

Se comparan ambas prediciones utilizando el indicador:

\[rmse = \sqrt{\frac{\sum_{i=i}^{N}(\widehat{y_{i}} - y_{i})^2}{N}}\]

rmse(M1[,1],M1[,6])
## [1] 6.742762
rmse(M2[,1],M2[,6])
## [1] 9.071189

Comportamiento del \(R_{2}\) :

# CONSISTENCIA R2
N = 10000


M = f_ProcesoGeneradorDatos(N,B0,B1,B2,B3,Sigma)

R = data.frame()


for(i in seq(10, N, by = 1)){

  P = cbind(f_RegresionLineal(M[1:i,]),"N"=i)

  R = rbind(R,P)
 
}


ggplot(data=R, aes(x=N, y=R2)) +
  geom_line(color="blue") + geom_hline(yintercept = 1, color = "red")

5. Predicción de un modelo que no cumple los supuestos

Se observa que los parámetros estimados no corresponden a los poblacionales, no se cumplen las propiedad asintóticas y las predicciones dentro y fuera de la muestra no se corresponden con los datos observados.

B0     = 10
B1     = 0.3
B2     = 0.4
B3     = 1.2
Sigma  = 8


# Datos de la muestra (entrenamiento)
M3 = f_ProcesoGeneradorDatos2(70,B0,B1,B2,B3,Sigma)

# Datos fuera de la muestra (test)
M4 = f_ProcesoGeneradorDatos2(30,B0,B1,B2,B3,Sigma)

# Regresion Lineal
mrl  = lm(Y~Z1+Z2+Z3,data = as.data.frame(M3))

summary(mrl)
## 
## Call:
## lm(formula = Y ~ Z1 + Z2 + Z3, data = as.data.frame(M3))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -94.20 -49.82 -22.92  16.63 382.69 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 77.82604   96.04657   0.810  0.42068   
## Z1           0.41922    0.58372   0.718  0.47517   
## Z2           0.28111    0.08277   3.396  0.00116 **
## Z3           1.27170    2.12055   0.600  0.55076   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 85.06 on 66 degrees of freedom
## Multiple R-squared:  0.1725, Adjusted R-squared:  0.1349 
## F-statistic: 4.585 on 3 and 66 DF,  p-value: 0.005619
confint(mrl)
##                    2.5 %      97.5 %
## (Intercept) -113.9371128 269.5892027
## Z1            -0.7462179   1.5846674
## Z2             0.1158486   0.4463655
## Z3            -2.9621182   5.5055260
# Prediccion dentro de la muestra
F = fitted(mrl,interval)
M3 = cbind(M3, Y_hat = F )


# Prediccion fuera de la muestra
P = predict(mrl, as.data.frame(M4[,3:5]),interval = 'confidence')
M4= cbind(M4, P )

rmse(M3[,1],M3[,6])
## [1] 82.59298
rmse(M4[,1],M4[,6])
## [1] 94.77064

Comportamiento del \(R_{2}\) :

# CONSISTENCIA R2
N = 10000


M = f_ProcesoGeneradorDatos2(N,B0,B1,B2,B3,Sigma)

R = data.frame()


for(i in seq(10, N, by = 1)){

  P = cbind(f_RegresionLineal(M[1:i,]),"N"=i)

  R = rbind(R,P)
 
}


ggplot(data=R, aes(x=N, y=R2)) +
  geom_line(color="blue") + geom_hline(yintercept = 1, color = "red")