Errores normales

• Se generan los valores de \(X\) que se utilizarán en todo el ejemplo.
• En este primer caso se consideran errores normales.

set.seed(1010)
x <- rnorm(176, 50, 15)
y <- -50 + 3.5*x + rnorm(176, 0, 25)

Exploración de los datos

• Gráfico de dispersión de los datos.

• La línea roja se agrega con el comando smooth.splines. Se utiliza para identificar tendencias en gráficos de dispersión.

• Se ajusta el modelo y se calculan los errores estandarizados con la expresión

\[ r_i = \frac{e_i}{\hat{\sigma}\sqrt{(1-h_i)}} \]

modelo <- lm(y ~ x)                   # Ajuste del modelo

res <- residuals(modelo)              # Residuos
s2h <- sum(res^2)/modelo$df.residual  # Varianza estimada

X <- cbind(1, x)                      # Matriz de diseño
H <- X%*%solve(t(X)%*%X)%*%t(X)       # Matriz sombrero

res_std <- res/(sqrt(1-diag(H))*sqrt(s2h)) # Residuos estandarizados

Métodos no estadísticos

• Histograma y densidad teórica

• Densidad estimada y densidad teórica

• La densidad se estima con density y se gráfica con plot

• Se validan las conclusiones anteriores con la simetría y curtosis

simetria <- function(x) mean((x-mean(x))^3)/sqrt(mean((x-mean(x))^2))^3
curtosis <- function(x) mean((x-mean(x))^4)/mean((x-mean(x))^2)^2 - 3

simetria(res_std)

## [1] 0.1715775

curtosis(res_std)

## [1] 0.09590063

• Gráfico cuantil-cuantil

• Para el gráfico se usa el comando qqnorm y para la recta qqline

Pruebas de normalidad

• Incluidas en el paquete stats

ks.test(res_std, 'pnorm') # Kolmogorov-Smirnov

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  res_std
## D = 0.044887, p-value = 0.8702
## alternative hypothesis: two-sided

shapiro.test(res_std)     # Shapiro-Wilk

## 
##  Shapiro-Wilk normality test
## 
## data:  res_std
## W = 0.99389, p-value = 0.6786

• Incluidas en el paquete nortest

library(nortest)
ad.test(res_std)          # Anderson-Darling

## 
##  Anderson-Darling normality test
## 
## data:  res_std
## A = 0.38506, p-value = 0.3893

cvm.test(res_std)         # Cramér-von Mises

## 
##  Cramer-von Mises normality test
## 
## data:  res_std
## W = 0.064456, p-value = 0.3288

lillie.test(res_std)      # Lilliefors

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  res_std
## D = 0.045327, p-value = 0.5068

pearson.test(res_std)     # Chi-cuadrada de Pearson

## 
##  Pearson chi-square normality test
## 
## data:  res_std
## P = 16.364, p-value = 0.23

Errores con colas pesadas

• Se generan los valores de \(Y\) con errores no normales

y2 <- -50 + 3.5*x + 11.6*rt(176, 5)

Exploración de los datos

• Gráfico de dispersión de los datos.

• Se ajusta el modelo y se calculan los errores estandarizados

modelo2 <- lm(y2 ~ x)                     # Ajuste del modelo

res2 <- residuals(modelo2)                # Residuos
s2h2 <- sum(res2^2)/modelo2$df.residual   # Varianza estimada

res_std2 <- res2/(sqrt(1-diag(H))*sqrt(s2h2)) # Residuos estandarizados

Métodos no estadísticos

• Histograma y densidad teórica

• Densidad estimada y densidad teórica

• Se validan las conclusiones anteriores con la simetría y curtosis

simetria(res_std2)

## [1] -0.02912022

curtosis(res_std2)

## [1] 1.803613

• Gráfico cuantil-cuantil

Pruebas de normalidad

• Incluidas en el paquete stats

ks.test(res_std2, 'pnorm') # Kolmogorov-Smirnov

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  res_std2
## D = 0.073817, p-value = 0.2929
## alternative hypothesis: two-sided

shapiro.test(res_std2)     # Shapiro-Wilk

## 
##  Shapiro-Wilk normality test
## 
## data:  res_std2
## W = 0.96749, p-value = 0.0003961

• Incluidas en el paquete nortest

ad.test(res_std2)          # Anderson-Darling

## 
##  Anderson-Darling normality test
## 
## data:  res_std2
## A = 1.5633, p-value = 0.0004885

cvm.test(res_std2)         # Cramér-von Mises

## 
##  Cramer-von Mises normality test
## 
## data:  res_std2
## W = 0.24705, p-value = 0.001381

lillie.test(res_std2)      # Lilliefors

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  res_std2
## D = 0.074732, p-value = 0.01798

pearson.test(res_std2)     # Chi-cuadrada de Pearson

## 
##  Pearson chi-square normality test
## 
## data:  res_std2
## P = 18.182, p-value = 0.1507

Errores uniformes

• Se generan los valores de \(Y\) con errores uniformes

y3 <- -50 + 3.5*x + runif(176, -25, 25)

Exploración de los datos

• Gráfico de dispersión de los datos.

• Se ajusta el modelo y se calculan los errores estandarizados

modelo3 <- lm(y3 ~ x)                     # Ajuste del modelo

res3 <- residuals(modelo3)                # Residuos
s2h3 <- sum(res3^2)/modelo3$df.residual   # Varianza estimada

res_std3 <- res3/(sqrt(1-diag(H))*sqrt(s2h3)) # Residuos estandarizados

Métodos no estadísticos

• Histograma y densidad teórica

• Densidad estimada y densidad teórica

• Se validan las conclusiones anteriores con la simetría y curtosis

simetria(res_std3)

## [1] -0.04002478

curtosis(res_std3)

## [1] -1.178003

• Gráfico cuantil-cuantil

Pruebas de normalidad

• Incluidas en el paquete stats

ks.test(res_std3, 'pnorm') # Kolmogorov-Smirnov

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  res_std3
## D = 0.088572, p-value = 0.1264
## alternative hypothesis: two-sided

shapiro.test(res_std3)     # Shapiro-Wilk

## 
##  Shapiro-Wilk normality test
## 
## data:  res_std3
## W = 0.95605, p-value = 2.616e-05

• Incluidas en el paquete nortest

ad.test(res_std3)          # Anderson-Darling

## 
##  Anderson-Darling normality test
## 
## data:  res_std3
## A = 1.9608, p-value = 5.135e-05

cvm.test(res_std3)         # Cramér-von Mises

## 
##  Cramer-von Mises normality test
## 
## data:  res_std3
## W = 0.26363, p-value = 0.0008718

lillie.test(res_std3)      # Lilliefors

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  res_std3
## D = 0.087766, p-value = 0.002128

pearson.test(res_std3)     # Chi-cuadrada de Pearson

## 
##  Pearson chi-square normality test
## 
## data:  res_std3
## P = 31.818, p-value = 0.002556

Errores asimétricos

• Se define la distribución doble exponencial asimétrica

rald <- function(n, theta, lambda){
  res <- c()
  for(i in 1:n)
  {
    aux <- runif(1) < theta
    lambda2 <- theta/(1-theta)*lambda
    res[i] <- aux*rexp(1, lambda) - (1-aux)*rexp(1, lambda2)
  }
  return(res)
}

• Se generan los valores de \(Y\) con errores asimétricos

y4 <- -50 + 3.5*x + rald(176, 0.8, 0.1)

Exploración de los datos

• Gráfico de dispersión de los datos.

• Se ajusta el modelo y se calculan los errores studentizados

modelo4 <- lm(y4 ~ x)                     # Ajuste del modelo

res4 <- residuals(modelo4)                # Residuos
s2h4 <- sum(res4^2)/modelo4$df.residual   # Varianza estimada

res_std4 <- res4/(sqrt(1-diag(H))*sqrt(s2h4)) # Residuos estandarizados

Métodos no estadísticos

• Histograma y densidad teórica

• Densidad estimada y densidad teórica

• Se validan las conclusiones anteriores con la simetría y curtosis

simetria(res_std4)

## [1] 1.176847

curtosis(res_std4)

## [1] 1.755852

• Gráfico cuantil-cuantil

Pruebas de normalidad

• Incluidas en el paquete stats

ks.test(res_std4, 'pnorm') # Kolmogorov-Smirnov

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  res_std4
## D = 0.11733, p-value = 0.01573
## alternative hypothesis: two-sided

shapiro.test(res_std4)     # Shapiro-Wilk

## 
##  Shapiro-Wilk normality test
## 
## data:  res_std4
## W = 0.92002, p-value = 3.084e-08

• Incluidas en el paquete nortest

ad.test(res_std4)          # Anderson-Darling

## 
##  Anderson-Darling normality test
## 
## data:  res_std4
## A = 4.0484, p-value = 4.116e-10

cvm.test(res_std4)         # Cramér-von Mises

## 
##  Cramer-von Mises normality test
## 
## data:  res_std4
## W = 0.69459, p-value = 6.727e-08

lillie.test(res_std4)      # Lilliefors

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  res_std4
## D = 0.11774, p-value = 3.117e-06

pearson.test(res_std4)     # Chi-cuadrada de Pearson

## 
##  Pearson chi-square normality test
## 
## data:  res_std4
## P = 39.818, p-value = 0.0001479

Solución al problema de no normalidad

• En estos ejemplos, no podemos corregir la no normalidad con transformaciones, porque la relación entre \(X\) y \(Y\) ya es lineal y transformar causaría problemas de no linealidad.

• Como la no normalidad afecta la confianza de los intervalos y la significancia de las pruebas, podemos hacer inferencias aproximadas basadas en bootstrap.

Intervalos de confianza

• Se ajusta el modelo \(Y = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}\) con las \(n\) observaciones originales. Se guardan los valores ajustados \(\hat{\mathbf{Y}}\) y los residuos \(\mathbf{e}\).

• Se selecciona una muestra aleatoria con reemplazo de tamaño \(n\) del los residuos \(\mathbf{e}\), la llamamos \(\mathbf{e}_{boot}\).

• Se calculan los valores de \(Y\) del bootstrap como \(\mathbf{Y}_{boot} = \hat{\mathbf{Y}} + \mathbf{e}_{boot}\).

• Se ajusta un modelo de regresión de \(\mathbf{Y}_{boot}\) contra las \(X\) originales. Se guarda \(\hat{\boldsymbol{\beta}}^{(1)}_{boot}\).

• Se repite lo anterior un número grande de veces \(M\) para obtener: \(\hat{\boldsymbol{\beta}}^{(1)}_{boot}, \hat{\boldsymbol{\beta}}^{(2)}_{boot}, \ldots, \hat{\boldsymbol{\beta}}^{(M)}_{boot}\). Los \(\hat{\beta}^{(j)}_{i, boot}\) constituyen una muestra bootstrap de la distribución de \(\hat{\boldsymbol{\beta}}\).

• El intervalo de confianza aproximados para \(\beta_i\) se calculan a partir de los cuantiles muestrales de \(\hat{\beta}^{(j)}_{i, boot}\), donde \(\hat{\beta}^{(j)}_{i, boot}\) denota la \(i\)-ésima entrada del vector \(\hat{\boldsymbol{\beta}}^{(j)}_{boot}\).

Pruebas de hipótesis

\[ H_0: \mathbf{A}\boldsymbol{\beta} = \mathbf{b} \qquad \text{vs.} \qquad H_1: \mathbf{A}\boldsymbol{\beta} \neq \mathbf{b} \]

• Ajustar el modelo original y guardar el valor del estadístico \(Q\) que se usaría para contrastar \(H_0\), lo llamamos \(Q_{obs}\). Por ejemplo, \(t\) para pruebas individuales, \(F\) para significancia conjunta o \(T^2\) para una hipótesis lineal general.

• Ajustar el modelo bajo \(H_0\). Por ejemplo, si \(H_0: \beta_1 = 0\), el modelo es \(Y = \beta_0 + \epsilon\). Se guardan los valores ajustados \(\hat{\mathbf{Y}}_0\) y los residuos \(\mathbf{e}_{0}\).

• Se selecciona una muestra aleatoria con reemplazo de tamaño \(n\) de los residuos \(\mathbf{e}_0\), la llamamos \(\mathbf{e}_{boot}\).

• Se calculan los valores de \(Y\) del bootstrap como \(\mathbf{Y}_{boot} = \hat{\mathbf{Y}}_0 + \mathbf{e}_{boot}\).

• Se ajusta un modelo de regresión \(\mathbf{Y}_{boot}\) contra las las \(X\) originales. Se guarda el estadístico \(t^{(j)}_i\).

• Se repite lo anterior un número grande de veces \(M\) para obtener: \(Q^{(1)}, Q^{(2)}, \ldots, Q^{(M)}\). Los \(Q^{(j)}\) constituyen una muestra bootstrap de la distribución de \(Q\), que se utilizará para contrastar \(H_0\).

• Para contrastar \(H_0: \beta_i = 0\) con una significancia \(\alpha\) se calcula el cuantil muestral \(t_i^{(1-\alpha/2)}\) y se compara con el estadístico obtenido en el modelo inicial, \(Q_{obs}\). Si \(\vert\,Q_{obs}\,\vert > Q^{(1-\alpha/2)}\), se rechaza \(H_0\).

Intervalos de confianza

Errores con colas pesadas

• Obtenemos la muestra bootstrap con el siguiente código (muestreamos simultáneamente para ambos parámetros)

beta2 <- list(beta0 = c(), beta1 = c()) 
M <- 2500              # Tamaño de muestra bootstrap
for (i in 1:M)
{
  errores_boot <- sample(res2, size = 176, replace = T)     # Muestra de errores estandarizados
  y_boot <- fitted(modelo2) + errores_boot                  # Cálculo de Y_boot
  modelo_boot <- lm(y_boot ~ x)                             # Ajuste del modelo
  beta2$beta0[i] <- coef(modelo_boot)[1]                    # Se guarda la estimación de beta0
  beta2$beta1[i] <- coef(modelo_boot)[2]                    # Se guarda la estimación de beta1
}

• Calculamos los intervalos de confianza 95% como los percentiles muestrales

lapply(beta2, quantile, probs = c(0.025, 0.975))

## $beta0
##      2.5%     97.5% 
## -53.22182 -36.13864 
## 
## $beta1
##     2.5%    97.5% 
## 3.220238 3.553822

• Calculamos los intervalos de confinaza asumiendo normalidad en los errores

confint(modelo2)

##                  2.5 %     97.5 %
## (Intercept) -53.076043 -36.149355
## x             3.224004   3.552581

Errores uniformes

• Obtenemos la muestra bootstrap con el siguiente código (muestreamos simultáneamente para ambos parámetros)

beta3 <- list(beta0 = c(), beta1 = c()) 
for (i in 1:M)
{
  errores_boot <- sample(res3, size = 176, replace = T)     # Muestra de errores estandarizados
  y_boot <- fitted(modelo3) + errores_boot                  # Cálculo de Y_boot
  modelo_boot <- lm(y_boot ~ x)                             # Ajuste del modelo
  beta3$beta0[i] <- coef(modelo_boot)[1]                    # Se guarda la estimación de beta0
  beta3$beta1[i] <- coef(modelo_boot)[2]                    # Se guarda la estimación de beta1
}

• Calculamos los intervalos de confianza 95% como los percentiles muestrales

lapply(beta3, quantile, probs = c(0.025, 0.975))

## $beta0
##      2.5%     97.5% 
## -54.21397 -39.65778 
## 
## $beta1
##     2.5%    97.5% 
## 3.298288 3.585009

• Calculamos los intervalos de confinaza asumiendo normalidad en los errores

confint(modelo3)

##                  2.5 %     97.5 %
## (Intercept) -54.198064 -39.221247
## x             3.295236   3.585963

Errores asimétricos

• Obtenemos la muestra bootstrap con el siguiente código (muestreamos simultáneamente para ambos parámetros)

beta4 <- list(beta0 = c(), beta1 = c()) 
for (i in 1:M)
{
  errores_boot <- sample(res4, size = 176, replace = T) # Muestra de errores estandarizados
  y_boot <- fitted(modelo4) + errores_boot                  # Cálculo de Y_boot
  modelo_boot <- lm(y_boot ~ x)                             # Ajuste del modelo
  beta4$beta0[i] <- coef(modelo_boot)[1]                    # Se guarda la estimación de beta0
  beta4$beta1[i] <- coef(modelo_boot)[2]                    # Se guarda la estimación de beta1
}

• Calculamos los intervalos de confianza 95% como los percentiles muestrales

lapply(beta4, quantile, probs = c(0.025, 0.975))

## $beta0
##      2.5%     97.5% 
## -46.00863 -36.47316 
## 
## $beta1
##     2.5%    97.5% 
## 3.383428 3.567276

• Calculamos los intervalos de confinaza asumiendo normalidad en los errores

confint(modelo4)

##                  2.5 %     97.5 %
## (Intercept) -46.132489 -36.550787
## x             3.383281   3.569279

Pruebas de hipótesis individuales: para \(\beta_0\)

• La hipótesis a contrastar es

\[ H_0: \beta_0 = 0 \qquad \text{vs.} \qquad H_1: \beta_0 \neq 0 \]

• El modelo bajo la \(H_0\) es

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

• El estadístico de prueba es el estadístico \(t\) calculado como

\[ t = \frac{\hat{\beta}_0}{\sqrt{\hat{V}(\hat{\beta}_0)}} \]

donde \(\hat{V}(\hat{\beta_0}) = \hat{\sigma}^2\left( \frac{1}{n} + \frac{\bar{x}^2}{S_{xx}} \right)\).

• Los estadísticos se pueden obtener a partir del summary de cada modelo ajustado.

Errores con colas pesadas

• Obtenermos y guardamos el estadístico \(t_{obs}\) del resumen del ajuste del modelo 2.

summary(modelo2)[[4]][1, 3]

## [1] -10.40389

• Se obtiene la muestra bootstrap con el siguiente código.

modelo2_0 <- lm(y2 ~ x - 1)     # Se ajusta el modelo bajo H0
res2_0 <- residuals(modelo2_0)  # Se calculan los residuos

t2_0 <- c()
for (i in 1:M)
{
  errores_boot <- sample(res2_0, size = 176, replace = T)   # Muestra de errores
  y_boot <- fitted(modelo2_0) + errores_boot                # Calculo de Y_boot
  fit <- lm(y_boot ~ x)                                     # Ajuste del modelo con la muestra bootstrap
  t2_0[i] <- summary(fit)[[4]][1, 3]                        # Se guarda el estadístico de prueba
}

• Se comparan los valores críticos

quantile(t2_0, probs = c(0.975))  # Cuantil bootstrap (no normalidad)

##    97.5% 
## 1.221038

qt(0.975, modelo2$df.residual)    # Cuantil teórico (bajo normalidad)

## [1] 1.973691

• En ambos casos la conclusión es rechazar \(H_0\).

Errores uniformes

• Obtenermos y guardamos el estadístico \(t_{obs}\) del resumen del ajuste del modelo 3.

summary(modelo3)[[4]][1, 3]

## [1] -12.31109

• Se obtiene la muestra bootstrap con el siguiente código.

modelo3_0 <- lm(y3 ~ x - 1)     # Se ajusta el modelo bajo H0
res3_0 <- residuals(modelo3_0)  # Se calculan los residuos

t3_0 <- c()
for (i in 1:M)
{
  errores_boot <- sample(res3_0, size = 176, replace = T)   # Muestra de errores
  y_boot <- fitted(modelo3_0) + errores_boot                # Calculo de Y_boot
  fit <- lm(y_boot ~ x)                                     # Ajuste del modelo con la muestra bootstrap
  t3_0[i] <- summary(fit)[[4]][1, 3]                        # Se guarda el estadístico de prueba
}

• Se comparan los valores críticos

quantile(t3_0, probs = c(0.975))  # Cuantil bootstrap (no normalidad)

##    97.5% 
## 1.035576

qt(0.975, modelo3$df.residual)    # Cuantil teórico (bajo normalidad)

## [1] 1.973691

• En ambos casos la conclusión es rechazar \(H_0\).

Errores asimétricos

• Obtenermos y guardamos el estadístico \(t_{obs}\) del resumen del ajuste del modelo 4.

summary(modelo4)[[4]][1, 3]

## [1] -17.03155

• Se obtiene la muestra bootstrap con el siguiente código.

modelo4_0 <- lm(y4 ~ x - 1)       # Se ajusta el modelo bajo H0
res4_0 <- residuals(modelo4_0)    # Se calculan los residuos

t4_0 <- c()
for (i in 1:M)
{
  errores_boot <- sample(res4_0, size = 176, replace = T)   # Muestra de errores
  y_boot <- fitted(modelo4_0) + errores_boot                # Calculo de Y_boot
  fit <- lm(y_boot ~ x)                                     # Ajuste del modelo con la muestra bootstrap
  t4_0[i] <- summary(fit)[[4]][1, 3]                        # Se guarda el estadístico de prueba
}

• Se comparan los valores críticos

quantile(t4_0, probs = c(0.975))  # Cuantil bootstrap (no normalidad)

##     97.5% 
## 0.9960963

qt(0.975, modelo4$df.residual)    # Cuantil teórico (bajo normalidad)

## [1] 1.973691

• En ambos casos la conclusión es rechazar \(H_0\).

Pruebas de hipótesis individuales: para \(\beta_1\)

• La hipótesis a contrastar es

\[ H_0: \beta_1 = 0 \qquad \text{vs.} \qquad H_1: \beta_1 \neq 0 \]

• El modelo bajo la \(H_0\) es

\[ Y = \beta_0 + \epsilon \]

• El estadístico de prueba es el estadístico \(t\) calculado como

\[ t = \frac{\hat{\beta}_1}{\sqrt{\hat{V}(\hat{\beta}_1)}} \]

donde \(\hat{V}(\hat{\beta_1}) = \frac{\hat{\sigma}^2}{S_{xx}}\).

• Obtenermos y guardamos el estadístico \(t_{obs}\) del resumen del ajuste del modelo 2.

summary(modelo2)[[4]][2, 3]

## [1] 40.70544

• Se obtiene la muestra bootstrap con el siguiente código.

modelo2_1 <- lm(y2 ~ 1)     # Se ajusta el modelo bajo H0
res2_1 <- residuals(modelo2_1)  # Se calculan los residuos

t2_1 <- c()
for (i in 1:M)
{
  errores_boot <- sample(res2_1, size = 176, replace = T)   # Muestra de errores
  y_boot <- fitted(modelo2_1) + errores_boot                # Calculo de Y_boot
  fit <- lm(y_boot ~ x)                                     # Ajuste del modelo con la muestra bootstrap
  t2_1[i] <- summary(fit)[[4]][2, 3]                        # Se guarda el estadístico de prueba
}

• Se comparan los valores críticos

quantile(t2_1, probs = c(0.975))  # Cuantil bootstrap (no normalidad)

##    97.5% 
## 1.934417

qt(0.975, modelo2$df.residual)    # Cuantil teórico (bajo normalidad)

## [1] 1.973691

• En ambos casos la conclusión es rechazar \(H_0\).

Errores uniformes

• Obtenermos y guardamos el estadístico \(t_{obs}\) del resumen del ajuste del modelo 3.

summary(modelo3)[[4]][2, 3]

## [1] 46.7152

• Se obtiene la muestra bootstrap con el siguiente código.

modelo3_1 <- lm(y3 ~ 1)     # Se ajusta el modelo bajo H0
res3_1 <- residuals(modelo3_1)  # Se calculan los residuos

t3_1 <- c()
for (i in 1:M)
{
  errores_boot <- sample(res3_1, size = 176, replace = T)   # Muestra de errores
  y_boot <- fitted(modelo3_1) + errores_boot                # Calculo de Y_boot
  fit <- lm(y_boot ~ x)                                     # Ajuste del modelo con la muestra bootstrap
  t3_1[i] <- summary(fit)[[4]][2, 3]                        # Se guarda el estadístico de prueba
}

• Se comparan los valores críticos

quantile(t3_1, probs = c(0.975))  # Cuantil bootstrap (no normalidad)

##    97.5% 
## 1.952551

qt(0.975, modelo3$df.residual)    # Cuantil teórico (bajo normalidad)

## [1] 1.973691

• En ambos casos la conclusión es rechazar \(H_0\).

Errores asimétricos

• Obtenermos y guardamos el estadístico \(t_{obs}\) del resumen del ajuste del modelo 4.

summary(modelo4)[[4]][2, 3]

## [1] 73.77609

• Se obtiene la muestra bootstrap con el siguiente código.

modelo4_1 <- lm(y4 ~ 1)           # Se ajusta el modelo bajo H0
res4_1 <- residuals(modelo4_1)    # Se calculan los residuos

t4_1 <- c()
for (i in 1:M)
{
  errores_boot <- sample(res4_1, size = 176, replace = T)   # Muestra de errores
  y_boot <- fitted(modelo4_1) + errores_boot                # Calculo de Y_boot
  fit <- lm(y_boot ~ x)                                     # Ajuste del modelo con la muestra bootstrap
  t4_1[i] <- summary(fit)[[4]][2, 3]                        # Se guarda el estadístico de prueba
}

• Se comparan los valores críticos

quantile(t4_1, probs = c(0.975))  # Cuantil bootstrap (no normalidad)

##    97.5% 
## 2.009705

qt(0.975, modelo4$df.residual)    # Cuantil teórico (bajo normalidad)

## [1] 1.973691

• En ambos casos la conclusión es rechazar \(H_0\).

Significancia conjunta

• La hipótesis a contrastar es

\[ H_0: \beta_1 = \ldots = \beta_p = 0 \qquad \text{vs.} \qquad H_1: \text{ algún }\beta_i \neq 0 \]

• El modelo bajo la \(H_0\) es

\[ Y = \beta_0 + \epsilon \]

En nuestro ejemplo, este modelo coincide con el modelo para contrastar \(H_0:\beta_1 = 0\).

• El estadístico de prueba es el estadístico \(F\) calculado como

\[ F = \frac{SCR/p}{SCE/(n-p-1)} \]

donde \(SCR\) es la suma de cuadrados de regresión \(\sum_{i=1}^n(\hat{y}_i-\bar{y}_n)^2\) y \(SCE\) es la suma de cuadrados del error \(\sum_{i=1}^n (y_i - \hat{y}_n)^2\).

Errores con colas pesadas

• Obtenermos y guardamos el estadístico \(F_{obs}\) del resumen del ajuste del modelo 2.

summary(modelo2)[[10]][1]

##    value 
## 1656.933

• Se obtiene la muestra bootstrap con el siguiente código.

F2 <- c()
for (i in 1:M)
{
  errores_boot <- sample(res2_1, size = 176, replace = T)   # Muestra de errores
  y_boot <- fitted(modelo2_1) + errores_boot                # Calculo de Y_boot
  fit <- lm(y_boot ~ x)                                     # Ajuste del modelo con la muestra bootstrap
  F2[i] <- summary(fit)[[10]][1]                            # Se guarda el estadístico de prueba
}

• Se comparan los valores críticos

quantile(F2, probs = c(0.95))       # Cuantil bootstrap (no normalidad)

##      95% 
## 4.125201

qf(0.95, 1, modelo2$df.residual)    # Cuantil teórico (bajo normalidad)

## [1] 3.895458

• En ambos casos la conclusión es rechazar \(H_0\).

Errores uniformes

• Obtenermos y guardamos el estadístico \(F_{obs}\) del resumen del ajuste del modelo 2.

summary(modelo3)[[10]][1]

##   value 
## 2182.31

• Se obtiene la muestra bootstrap con el siguiente código.

F3 <- c()
for (i in 1:M)
{
  errores_boot <- sample(res3_1, size = 176, replace = T)   # Muestra de errores
  y_boot <- fitted(modelo3_1) + errores_boot                # Calculo de Y_boot
  fit <- lm(y_boot ~ x)                                     # Ajuste del modelo con la muestra bootstrap
  F3[i] <- summary(fit)[[10]][1]                            # Se guarda el estadístico de prueba
}

• Se comparan los valores críticos

quantile(F3, probs = c(0.95))       # Cuantil bootstrap (no normalidad)

##     95% 
## 3.93877

qf(0.95, 1, modelo3$df.residual)    # Cuantil teórico (bajo normalidad)

## [1] 3.895458

• En ambos casos la conclusión es rechazar \(H_0\).

Errores asimétricos

• Obtenermos y guardamos el estadístico \(F_{obs}\) del resumen del ajuste del modelo 4.

summary(modelo4)[[10]][1]

##    value 
## 5442.911

• Se obtiene la muestra bootstrap con el siguiente código.

F4 <- c()
for (i in 1:M)
{
  errores_boot <- sample(res4_1, size = 176, replace = T)   # Muestra de errores
  y_boot <- fitted(modelo4_1) + errores_boot                # Calculo de Y_boot
  fit <- lm(y_boot ~ x)                                     # Ajuste del modelo con la muestra bootstrap
  F4[i] <- summary(fit)[[10]][1]                            # Se guarda el estadístico de prueba
}

• Se comparan los valores críticos

quantile(F4, probs = c(0.95))       # Cuantil bootstrap (no normalidad)

##      95% 
## 3.809689

qf(0.95, 1, modelo4$df.residual)    # Cuantil teórico (bajo normalidad)

## [1] 3.895458

• En ambos casos la conclusión es rechazar \(H_0\).