Tabajo Práctico N°2
Estadística Aplicada II
Ignacio Santamaria
101613
Caso de Aplicación
Una empresa de retail, desea conocer el impacto en las ventas los catálogos de productos enviados a los hogares de sus clientes. Para eso diseña 4 estrategias: realizar envíos cada 1, 2, 3 y 4 meses y recolecta la información durante dos años, incluyendo otras variables que considera que pueden ser relevantes:
Parte I:
Teniendo en cuenta que su objetivo es predecir el Monto anual gastado, y utilizando como variable explicativa únicamente el Salario:
library(readxl)## Warning: package 'readxl' was built under R version 4.0.3datos_parte1 <- read_excel("C:/Users/Ignac/Desktop/Estadistica Aplicada II/TPs/TP2/Parte I.xlsx") #importo datos del excel
head(datos_parte1)a) Estime el modelo lineal. Construya el diagrama de dispersión y el de residuos. Interprete e indique que inconvenientes observa en este modelo.
Panteo un modelo de regresión lineal simple (RLS), utilizando las siguientes variables:
\(Y_i\): Monto gastado en compras, por el i-ésimo cliente en el último año (montoAnual)[$]. Variable aleatoria, de respuesta/explicada.
\(X_i\): Salario del i-ésimo cliente (salario)[$]. Variable considerada determinística, explicativa.
Modelo:
\(Y_i = \beta_0 +\beta_1*X_i+\tilde{\epsilon}_i\) recta de regresión poblacional
Siendo \(\tilde{\epsilon}\) la perturbación aleatoria.
\(Y_i = b_0 + b_1*X_i\) recta de regresión estimada
Estimo los coeficientes de la recta, b0 y b1, mediante el método de mínimos cuadrados, utilizando las siguientes sumas de cuadrados:
\(SC_x=\displaystyle{\sum^{n}_{i=1} (X_i-\bar{X})^2}=\displaystyle{\sum^{n}_{i=1} X_i^2-n*\bar{X}^2}\)
\(SC_{xy}=\displaystyle{\sum^{n}_{i=1} (X_i-\bar{X})*(Y_i-\bar{Y})}=\displaystyle{\sum^{n}_{i=1} (X_i*Y_i)-n*\bar{X}*\bar{Y}}\)
=> \(b_1=\frac{SC_{xy}}{SC_x}\) y \(b_0=\bar{Y}-b_1*\bar{X}\)
options(scipen = 999)
salario <- datos_parte1$salario #observaciones de Y
montoAnual <- datos_parte1$montoAnual #observaciones de X
sc_x <- sum((salario-mean(salario))^{2}) #suma de cuadrados de las x
sc_xy <- sum((salario-mean(salario))*(montoAnual-mean(montoAnual))) #suma de cuadrados de las y
b1 <- sc_xy/sc_x #estimación de la pendiente de la recta
b0 <- mean(montoAnual)-b1*mean(salario) #estomación de la ordenada al origen del al recta
print(c(b1=b1, b0=b0))## b1 b0
## 2.880729 -42542.095535\(Y_i = -42542.0955[\$monto] + 2.8807[\frac{\$salario}{\$monto}]*X_i\) Estimación del modelo lineal de regresión simple.
Para este modelo de RLS con los datos dados, se puede interpretar el coeficiente de regresión estimado b1=2.8807 (pendiente de la recta) como que por cada peso (unidad monetaria $) que aumente el salario del cliente, entonces el monto gastado por este mismo aumentará 2.8807 pesos. En cuanto a b0=-42542.0955 (ordenada al origen), para este caso, este valor no tiene una interpretación real más allá de la matemática. Ya que no tendría sentido considerar un cliente con salario nulo y menos aún pensar que tendría un monto anual gastado negativo.
Grafico el diagrama de dispersión de los datos, la recta con los coeficientes estimados y marco la media muestral (“centro de masa”) de la variable salario con una línea a rayas:
library(ggplot2)
ggplot(datos_parte1, aes(x=salario, y=montoAnual)) +
geom_point(color="blue") + # armo el gráfico de dispersión de los datos
labs(x="Salario del cliente", y="Monto anual gastado por el cliente",
title="Diagrama de dispersión de datos observados",
subtitle="con recta de regresión estimada") +
geom_smooth(method=lm, color="green", se=FALSE) + # grafico la recta de regresión estimada
geom_label(label="8", x=100000, y=580000, label.size=0, fill="lightblue") +
geom_label(label="84", x=118000, y=630000, label.size=0, fill="lightblue") +
geom_label(label="88", x=128000, y=600000, label.size=0, fill="lightblue") +
geom_vline(xintercept=mean(salario), linetype="dashed") # promedio muestral de las X (línea vertical)
Calculo los predichos \(\hat{Y_i}\) y los residuos \(e_i\)con las siguientes expresiones:
\(\hat{Y_i}=b_0 + b_1*X_i\) y \(e_i=Y_i-\hat{Y_i}\)
y_pico <- -42542.095535 + 2.880729*salario # calculo los predichos
e <- montoAnual - y_pico # calculo los residuosArmo un gráfico de residuos vs predichos:
ggplot(data = data.frame("predichos"=y_pico, "residuos"=e), aes(x=y_pico, y=e)) +
geom_point(color="orange") +
labs(x="Monto anual gastado estimado (predichos)", y="Residuos",
title="Residuos VS Predichos") +
geom_hline(yintercept=0, linetype="dashed", color="red")
Análisis de Especificación a partir de los diagramas:
En primera instancia se observa en el diagrama de dispersión las observaciones de la variable X (salario del cliente) y de la variable Y (monto anual gastado por el cliente), junto con el gráfico de la recta estimada por cuadrados mínimos. A simple vista, un primer llamado de atención podria ser que hay varios puntos (para los cuales marqué los números de las respectivas observaciones) que están muy alejados. Tanto de la recta, por lo que su discrepancia debe ser alta (residuos de valor posiblemente alto), como así tambien están alejados del “centro de masa” de las Y (promedio muestral), por lo que su leverage tambien será posiblemente alto. Esto podría llevar a pensar en que serán posibles puntos influyentes y por ende posibles outliers.
Sin embargo, observando la nube de puntos completa se ve que a mayores valores de X e Y los puntos tienden a formar una curva que crece exponencialmente, y además esta nube de puntos se expande cada vez más. Esta debe ser la posible causa de que se observe una posible alta discrepancia y un posible alto leverage en los puntos antes mencionados. Entonces, el verdadero inconveniente que se observa es que la relación entre las variables X e Y NO es lineal, sino que es mas bien exponencial creciente o polinómica de un término.
En segundo lugar se observa el gráfico de residuos VS predichos. Ya que los residuos son estimaciones puntuales de las perturbaciones aleatorias, entonces en este gráfico debería observarse que se cumplan los supuestos de Ausenia de Vicio (\(E[\tilde{\epsilon_i}]=0\)) y de Homocedasticidad (\(D^{2}[\tilde{\epsilon_i}]=\sigma^{2}\), varianza constante).
Para esta nube de puntos se puede ver que (similar al caso del diagrama de dispersión) la nube se expande abruptamente a medida que aumentan los valores de los predichos. Esto indica que claramente NO se cumple el supuesto de Homocedasticidad de las perturbaciones ya que aumenta mucho la variabilidad de los residuos a medida que aumenta el valor de los predichos, y por lo tanto la varianza no va a ser constante.
Por otra parte, para la parte de la nube de puntos de valores de predichos bajos (cercanos al predicho de valor cero) se puede ver claramente que el conjunto de puntos esta absolutamente todo por encima de la línea de valor nulo de los residuos. Esto indica que NO se cumple el supuesto de Ausencia de Vicio porque los residuos NO varian aleatoriamente alrededor de cero, y por ende NO se puede decir que la media de las perturbaciones sea nula.
En conclusión, la especificación del modelo presentado NO es correcta, tal cual se ve evidenciado en los residuos.
b) Identifique alguna transformación que le permita resolver el/los inconvenientes/s anterior/es.
Para resolver el inconveniente en la especificación de este modelo de regresión, y dado que se observó, en el diagrama de dispersión, que la relación entre el monto anual gastado por el cliente y el salario de este mismo es posiblemente exponencial creciente o polinómica de un término, se propone plantear Y en función de X de esta primera forma.
Modelo propuesto (no lineal):
\(Y_i=\theta_0*e^{{\theta_1}*X_{i}}*\epsilon'_i\)
Entonces, como transformación linealizante puedo aplicar logaritmo natural, obteniendo una transformación semi-logaritmica, tal cual sugeriría la Regla de Mosteller y Tukey. Por lo que se transformará solo la variable Y (si hubiese tomado una función polinómica en vez de exponencial como modelo, hubiese obtenido una transformación logarítmica, habiese transformado ambas variables).
\(ln(Y_i)=ln(\theta_0*e^{{\theta_1}*X_{i}}*\epsilon'_i)\)
\(ln(Y_i)=ln(\theta_0)+\theta_1*X_i+ln(\epsilon'_i)\)
=> \(\beta_0=ln(\theta_0)\) , \(\beta_1=\theta_1\) y \(\tilde{\epsilon}_i=ln(\epsilon'_i)\)
Finalmente se llega al siguiente modelo linealizado de regresión lineal simple:
\(ln(Y_i)=\beta_0+\beta_1*X_i+\tilde{\epsilon}_i\)
Ahora transformo las muestras dadas para las observaciones de X e Y. Y como se hizo en el inciso a), calculo los coeficientes de regresión y grafico los diagramas de dispersión, la recta estimada y los residuos vs predichos:
datos_semi_transformados_parte1 <- data.frame("salario"=salario, "ln_montoAnual"=log(montoAnual)) # datos transformados
salario <- datos_semi_transformados_parte1$salario
ln_montoAnual <- datos_semi_transformados_parte1$ln_montoAnual
sc_ln_x <- sum((salario-mean(salario))^{2}) #suma de cuadrados de las X
sc_ln_xy <- sum((salario-mean(salario))*(ln_montoAnual-mean(ln_montoAnual))) #suma de cuadrados de las ln(Y)
b1_ln <- sc_ln_xy/sc_ln_x #estimación de la pendiente de la nueva recta
b0_ln <- mean(ln_montoAnual)-b1_ln*mean(salario) #estimación de la ordenada al origen del la nueva recta
print(c(b1=b1_ln, b0=b0_ln))## b1 b0
## 0.00002190898 10.03612280857Entonces, con la transformación semi-logaritmica, la estimación de la recta de regresión queda:
\(Y_i = 10.03612280857[\$monto] + 0.00002190898[\frac{\$salario}{\$monto}]*X_i\) Estimación del modelo lineal de regresión simple habiendo aplicado la transformación semi-logarítmica.
ggplot(datos_semi_transformados_parte1, aes(x=salario, y=ln_montoAnual)) +
geom_point() + # armo el gráfico de dispersión de los datos
labs(x="Salario del cliente", y="ln(Monto anual gastado por el cliente)",
title="Diagrama de dispersión de datos observados",
subtitle="con recta de regresión estimada") +
geom_smooth(method=lm, color="darkgreen", se=FALSE) + # grafico la recta de regresión estimada
geom_vline(xintercept = mean(salario), linetype="dashed") # promedio muestral de las X (línea vertical)## `geom_smooth()` using formula 'y ~ x'
y_pico_ln <- 10.03612280857 + 0.00002190898*salario # calculo los predichos
e_ln <- ln_montoAnual - y_pico_ln # calculo los residuos
ggplot(data = data.frame("predichos"=y_pico_ln, "residuos"=e_ln), aes(x=y_pico_ln, y=e_ln)) +
geom_point(color="red") +
labs(x="ln(Monto anual gastado estimado (predichos))", y="Residuos",
title="Residuos VS Predichos") +
geom_hline(yintercept=0, linetype="dashed")
Ahora se ve graficamente que gracias a la transformación aplicada se resuelven los inconvenientes anteriores. Ya que en el diagrama de dispersión la nube de puntos tiene una forma mucho mas lineal y NO presenta características curvilíneas. Mientras que observando el diagrama de Residuos VS Predichos, se puede decir que la variabilidad de los puntos se mantiene constante, pudiendose meter estos puntos dentro de una banda de variaciones uniforme y por ende cumpliendo el supuesto de Homocedasticidad. Y finalmente, en este mismo gráfico se ve que los puntos varían aleatoriamente respecto del valor nulo de los residuos y por lo tanto se cumple el supuesto de Auscencia de Vicio.
Parte II
Para cada uno de los siguientes casos, realice una codificación adecuada al tipo de variable, estime los modelos de regresión y responda las siguientes preguntas:
- II.a Variable Dependiente: ln(montoAnual) Independientes: ln(salario) y region
- II.b Variable Dependiente: ln(montoAnual) Independientes: ln(salario) y frecuencia
Para cada uno de estos modelos conteste:
- Interprete en términos del problema el significado de los coeficientes de regresión. Realice esta interpretación sean o no estos coeficientes significativos.
- Indique que coeficientes son significativos e indique que modelo final elegiría si solo contara con estas variables.
Observación: en esta parte se le aplicó logaritmo natural tanto a la variable montoAnual como a la variable salario.
a) Variable Dependiente: ln(montoAnual) Independientes: ln(salario) y region.
library(readxl)
Parte_2_A <- read_excel("C:/Users/Ignac/Desktop/Estadistica Aplicada II/TPs/TP2/Parte II a.xlsx")
head(Parte_2_A) # visualizo las primeras 6 filasDefino las siguientes variables:
\(Y_i\):(ln(montoAnual)) Logaritmo natural del monto anual gastado por i-ésimo cliente [$]. Variable aleatoria, de respuesta/explicada.
\(X_{1i}\):(ln(salario)) Logaritmo natural del salario del i-ésimo cliente [$]. Variable considerada determinística, explicativa.
\(X_{2i}\):(region) región de la localidad donde reside el i-ésimo cliente [CUALITATIVA]. Variable considerada determinística, explicativa.
Aplico logaritmo natural a las observaciones de las variables salario y montoAnual:
library(dplyr)
Parte_2_A_ln <- transmute(Parte_2_A, ln_salario = log(salario), region = region, ln_montoAnual = log(montoAnual))
head(Parte_2_A_ln) # visualizo las primeras 6 filas de los datos habiendo aplicado el logaritmo naturalArmo un gráfico de dispersión de los datos, separando por colores y formas los puntos correspondientes a cada categoria de la variable region para un mejor entendimiento de los mismos:
library(ggplot2)
ggplot(Parte_2_A_ln, aes(x = ln_salario, y = ln_montoAnual, shape = region, color = region)) +
geom_point() +
labs(title="Diagrama de dispersión", subtitle="con la variable categórica Región",
x="ln(Salario del cliente)", y="ln(Monto anual gastado por el cliente)", shape="Región", colour="Región")
Codificación de la variable categórica, especificación del modelo de RLM e interpretación de los coeficientes de regresión:
Dado que se tiene una variable cualitativa, region, y teniendo C = 4 categorias posibles para esta variable (Norte, Sur, Este, Oeste), voy a plantear (C-1) = 3 variables dummys (/binarias/indicadoras). Tomando como base/referancia la categoria Sur, se definen las siguiente varaibles D2, D3 y D4:
dummys_A <- matrix(data = c(0,0,0,1,0,0,0,1,0,0,0,1), byrow=TRUE, ncol=3, nrow=4,
dimnames=list(c("Sur", "Norte", "Este", "Oeste"), c("D2", "D3", "D4")))
print(dummys_A)## D2 D3 D4
## Sur 0 0 0
## Norte 1 0 0
## Este 0 1 0
## Oeste 0 0 1Representado en el dataframe:
region <- Parte_2_A_ln$region
D2 <- c()
for (i in 1:length(region)) {
ifelse(region[i] == "Norte", # la D2 corresponde a la categoría Norte
D2 <- c(D2, 1), D2 <- c(D2, 0))
}
D3 <- c()
for (i in 1:length(region)) {
ifelse(region[i] == "Este", # la D3 corresponde a la categoría Este
D3 <- c(D3, 1), D3 <- c(D3, 0))
}
D4 <- c()
for (i in 1:length(region)) {
ifelse(region[i] == "Oeste", # la D4 corresponde a la categoría Oeste
D4 <- c(D4, 1), D4 <- c(D4, 0))
}
Parte_2_A_ln_dummys <- data.frame("ln_salario"=Parte_2_A_ln$ln_salario,
"region"=region,"ln_montoAnual"=Parte_2_A_ln$ln_montoAnual,
"D2"=D2, "D3"=D3, "D4"=D4)
head(Parte_2_A_ln_dummys)Saco la variable categórica region porque ya no es necesaria:
ln_salario <- log(salario)
Parte_2_A_ln_dummys <- data.frame("ln_salario"=ln_salario, "ln_montoAnual"=ln_montoAnual,
"D2"=D2, "D3"=D3, "D4"=D4)Por lo tanto planteo el siguiente modelo de regresión lineal múltiple (RLM) usando las dummys:
\(Y_i=\beta_0+\beta_1*X_{1i}+\beta_2*D_2+\beta_3*D_3+\beta_4*D_4+\tilde{\epsilon_i}\)
Especificando el modelo de esta forma se logra modelar el efecto sobre la ordenada al origen (efecto aditivo) que genera la variable categórica region. Es decir, estos coeficientes se pueden interpretar como que independientemente del valor de la variable ln(salario) (es decir, independientemente del salario del cliente) (\(X_{1i}\)), la media de la variable ln(montoAnual) \(Y_i\) se va a incrementar en \(\beta_k\) (k=[2;4]) con el valor de la categoría correspondiente a \(D_k=1\), es decir que la ordenada al origen quedaría \(\beta_0+\beta_k\).
Aunque habría que análizar en más profundidad este caso de aplicación, ya que se desconoce de que tipo de empresa de retail se está hablando (de comida, indumentaria, etc.), se podría suponer que según la región en la que se encuentre el cliente, su deseo o necesidad de consumo de los productos de retail puede variar por las condiciones de la región, como por ejemplo condiciones sociales y culturales de moda o costumbres de consumo. Entonces, también puedo considerar en el modelo un posible efecto sobre la pendiente (interacción entre las dummys y la variable cuantitativa X1) que genera la variable categórica region, porque puede que cada dummy vaya a influir diferente sobre la variable de respuesta ln(montoAnual) \(Y_i\) dependiendo del valor de la variable cuantitativa \(X_{1i}\). Por lo que para representar esta interacción en el modelo debo especificarlo del siguiente modo:
\(Y_i=\beta_0+\beta_1*X_{1i}+\beta_2*D_2+\beta_3*D_3+\beta_4*D_4+\beta_5*(D_2*X_{1i})+\beta_6*(D_3*X_{1i})+\beta_7*(D_4*X_{1i})+\tilde{\epsilon_i}\)
Entonces, estos nuevos coeficientes se pueden interpretar como: que el aumento (pendiente \(\beta_1\)) de la variable de respuesta ln(montoAnual) respecto de cada valor que aumente ln(salario), se va a incrementar en \(\beta_m\) (m=[5;7]) con la categoría correspondiente a \(D_{(m-3)}=1\), es decir, que la pendiente del término de la variable ln(salario) quedaría \((\beta_1 + \beta_m)*X_{1i}\).
Al anterior dataframe ahora debo incluirle los datos de las variables \((D_k*X_{1i})\) (k=[2;4]). Quedando los datos del siguiente modo:
library(dplyr)
Parte_2_A_ln_dummys <- mutate(Parte_2_A_ln_dummys,
D2_mult_X = D2 * ln_salario, # Agrego las nuevas columnas de interacciones
D3_mult_X = D3 * ln_salario,
D4_mult_X = D4 * ln_salario)
head(Parte_2_A_ln_dummys)Con los datos anteriores se arma la matriz de datos para este modelo.
En cuanto al coeficiente \(\beta_1\), correspondiente a la variable explicativa cuantitativa ln(salario) (\(X_{1i}\)), su interpretación es que la media de la variable de respuesta ln(montoAnual) va a aumentar en \(\beta_1\) por cada unidad que aumente \(X_{1j}\). Esto es siempre manteniendoce CONSTANTES el resto de las variables explicativas (Ceteris Paribus).
Para que se entiendan mejor estas interpretaciones, expreso cómo quedaría la media de ln(montoAnual) para cada categoría posible de la variable region:
Si la categoría es Sur: \(E[Y|D_2=0,D_3=0,D_4=0]=\beta_0+\beta_1*X_1\)
Si la categoría es Norte: \(E[Y|D_2=1,D_3=0,D_4=0]=(\beta_0+\beta_2)+(\beta_1+\beta_5)*X_1\)
Si la categoría es Este: \(E[Y|D_2=0,D_3=1,D_4=0]=(\beta_0+\beta_3)+(\beta_1+\beta_6)*X_1\)
Si la categoría es Oeste: \(E[Y|D_2=0,D_3=0,D_4=1]=(\beta_0+\beta_4)+(\beta_1+\beta_7)*X_1\)
Ahora que tengo mi modelo RLM armado con los datos necesarios, calculo los coeficientes de regresión con los cálculos matriciales:
# Armo el VECTOR DE OBSERVACIONES DE Y
ln_montoAnual <- Parte_2_A_ln_dummys$ln_montoAnual
ln_montoAnual <- matrix(data = ln_montoAnual, nrow = length(ln_montoAnual))
# Armo la MATRIZ DE DATOS de dimensión (n x p)
n <- 100 # cantidad de obs de cada varaibles
p <- 8 # cantidad de parámetros a estimar
unos <- c()
for (i in 1:n) {
unos <- c(unos, 1) # creo la columna de unos
}
ln_salario <- Parte_2_A_ln_dummys$ln_salario # columna de ln(salario)
D2_mult_X <- Parte_2_A_ln_dummys$D2_mult_X # columna de D2*X
D3_mult_X <- Parte_2_A_ln_dummys$D3_mult_X # columna de D3*X
D4_mult_X <- Parte_2_A_ln_dummys$D4_mult_X # columna de D4*X
matriz_datos <- matrix(data = c(unos, ln_salario, D2, D3, D4, D2_mult_X, D3_mult_X, D4_mult_X), nrow=n, ncol=p)
head(matriz_datos) # muestro las primeras 6 filas de la matriz de datos## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 9.532424 0 0 0 0.00000 0.00000 0
## [2,] 1 11.247657 0 1 0 0.00000 11.24766 0
## [3,] 1 10.515967 0 0 0 0.00000 0.00000 0
## [4,] 1 11.415313 0 0 0 0.00000 0.00000 0
## [5,] 1 11.513925 0 1 0 0.00000 11.51392 0
## [6,] 1 10.599132 1 0 0 10.59913 0.00000 0Ahora armo la matriz \((X^{T}*X)^{-1}\), a la cual la voy a llamar “matriz normal” porque proviene de la ecuación normal matricial que resuelve este sistema:
matriz_normal <- solve(t(matriz_datos) %*% matriz_datos)Calculo el Vector de Coeficientes Estimado \(\beta=(X^{T}*X)^{-1}*X^T*Y\).
vect_coef_estimado <- matriz_normal %*% t(matriz_datos) %*% ln_montoAnual
rownames(vect_coef_estimado) <- c("Beta0", "Beta1", "Beta2", "Beta3", "Beta4", "Beta5", "Beta6", "Beta7")
colnames(vect_coef_estimado) <- c("Estimación Puntual")
vect_coef_estimado## Estimación Puntual
## Beta0 1.2504676
## Beta1 0.9221962
## Beta2 -0.9810033
## Beta3 -1.4445821
## Beta4 -4.2247655
## Beta5 0.1097512
## Beta6 0.1386216
## Beta7 0.4134926Se puede observar en el siguiente gráfico de dispersión, a las rectas estimadas para cada categoría de la variable region. Es decir, estas rectas representan la estimación de la media de la variable de respuesta ln(montoAnual) dada una categoria de la variable region, osea que son las estimaciones de las rectas que expresé anteriormente:
ggplot(Parte_2_A_ln, aes(x = ln_salario, y = ln_montoAnual, shape = region, color = region)) +
geom_point() +
labs(title="Diagrama de dispersión", subtitle="con rectas estimadas para cada categoría de región",
x="ln(Salario del cliente)", y="ln(Monto anual gastado por el cliente)", shape="Región", colour="Región") +
geom_smooth(method = lm, se = FALSE)
Tests de significación para los coeficientes y elección del modelo
Se realizará inferencia sobre los parámetros de regresión, utilizando el siguiente estadístico:
\(t_{n-p}=\frac{b_j-\beta_j}{\hat{D}(b_j)}\) con \(\hat{D}(b_j)=S*\sqrt{C_{ii}}\)
Siendo \(C_{ii}\) los elementos de la diagonal de la llamada “matriz normal”.
Test de hipótesis UNIlateral para el coeficiente \(\beta_1\):
\(H_0)\beta_1 \leq 0\) ; \(H_a)\beta_1 > 0\)
Dado que se puede suponer extra-estadísticamente que el ln(montoAnual) debería aumentar si aumenta ln(salario) entonces la hipótesis alternativa para este coeficiente será que este tenga que ser mayor que cero (pendiente positiva).
Test de hipótesis BIlateral para el resto de los 7 coeficientes estimados:
\(H_0)\beta_j=0\) ; \(H_a)\beta_j \neq 0\)
\((j \neq 1)\)
Entonces para calcular los valores P respectivos para cada uno de estos test, necesito calcular el desvío residual y el desvío de cada coeficiente estimado, y además voy a calcular un intervalo de confianza del 95% para cada coeficiente:
predichos <- matriz_datos %*% vect_coef_estimado # predichos para los datos observados
Q <- sum((ln_montoAnual-predichos)^{2}) # suma de cuadrado de los residuos
S <- sqrt(Q/(n-p)) # desvío residual
# Calculo el desvío de cada coeficiente estimado
C_ii <- diag(matriz_normal) # elementos de la diagonal de la "matriz normal"
desvio_coef_estimado <- S*sqrt(C_ii)
# Calculo los valores p (los calculo todos igual y despues modifico el de beta1)
vect_tobs <- (vect_coef_estimado-0)/desvio_coef_estimado # estadísticos t observados
valoresP <- 2*pt(abs(vect_tobs), df=n-p, lower.tail=FALSE) # valores P test BIlateral
valoresP[2] = valoresP[2]/2 # modifico el valor p para el beta 1 (test unilateral)
# Calculo los IC del 95% (NC=0.95=1-alfa=1-0.05)
LI_coef_estimado <- vect_coef_estimado-qt(1-0.05/2, df=n-p)*desvio_coef_estimado
LS_coef_estimado <- vect_coef_estimado+qt(1-0.05/2, df=n-p)*desvio_coef_estimado
# resultados de los test de significación
resultados <- c("---", "SIGNIFICATIVO", "NO significativo", "NO significativo", "**", "NO significativo", "NO significativo", "**")
matrix(data = c(vect_coef_estimado, desvio_coef_estimado, LI_coef_estimado, LS_coef_estimado, vect_tobs, valoresP, resultados),
nrow = 8, ncol = 7,
dimnames = list(c("Beta0", "Beta1", "Beta2", "Beta3", "Beta4", "Beta5", "Beta6", "Beta7"),
c("Estimación puntual", "Desvío estimado del coeficiente", "LI (IC95%)", "LS (IC95%)", "t observado", "Valor P", "Resultado")))## Estimación puntual Desvío estimado del coeficiente
## Beta0 "1.25046758818106" "1.91635448986071"
## Beta1 "0.922196226794961" "0.178663602647693"
## Beta2 "-0.981003280621252" "3.10690845115699"
## Beta3 "-1.44458205654254" "2.77805156314383"
## Beta4 "-4.224765517258" "2.77678884065729"
## Beta5 "0.109751177675853" "0.286301295992468"
## Beta6 "0.138621636820633" "0.256538769912494"
## Beta7 "0.413492605640821" "0.256533783032548"
## LI (IC95%) LS (IC95%) t observado
## Beta0 "-2.55557784255915" "5.05651301892128" "0.652524151871273"
## Beta1 "0.567354890239184" "1.27703756335074" "5.16163456422313"
## Beta2 "-7.15159164348398" "5.18958508224148" "-0.315749014186734"
## Beta3 "-6.9620322538872" "4.07286814080213" "-0.519998287903537"
## Beta4 "-9.73970783875005" "1.29017680423404" "-1.52145725141201"
## Beta5 "-0.458867908820162" "0.678370264171869" "0.383341532895963"
## Beta6 "-0.370886503870046" "0.648129777511313" "0.540353556961068"
## Beta7 "-0.0960056306758316" "0.922990841957474" "1.6118446496708"
## Valor P Resultado
## Beta0 "0.51569024972287" "---"
## Beta1 "0.000000702481075192944" "SIGNIFICATIVO"
## Beta2 "0.752908196419748" "NO significativo"
## Beta3 "0.604313781353778" "NO significativo"
## Beta4 "0.131574864732499" "**"
## Beta5 "0.702350931461651" "NO significativo"
## Beta6 "0.590258974809908" "NO significativo"
## Beta7 "0.110421256862023" "**"NOTA: el valor P correspondiente al \(\beta_1\) es para test UNIlateral, fue modificado dado \(\alpha^{\star}_{unilat}=\frac{\alpha^{\star}_{bilat}}{2}\). El resto de los valores P sin para tests BIlaterales.
\(Y_i=\beta_0+\beta_1*X_{1i}+\beta_2*D_2+\beta_3*D_3+\beta_4*D_4+\beta_5*(D_2*X_{1i})+\beta_6*(D_3*X_{1i})+\beta_7*(D_4*X_{1i})+\tilde{\epsilon_i}\)
Por lo tanto, si el test de significación para un coeficiente tiene un vapor P lo suficientemente bajo, entonces se rechazará la hipótesis nula del test con el respectivo valor P, ya que este será la probabilidad de riesgo de equivocarse al rechazar la hipótesis nula. Es decir, este valor P es la probabilidad de equivocarse al decir que la variable, a la que se corresponde el coeficiente, es significativa y que aporta información al modelo.
Entonces, se tiene que el \(\beta_1\) de la variable ln(salario) \(X_{1i}\) es significativo con solo un 0.00007025% de probabilidad de equivocarse. Otra forma de ver este resultado es diciendo que las variables ln(montoAnual) y ln(salario) presentan un cierto grado de asociación lineal, con solo 0.00007025% de probabilidad de equivocarse.
** Los coeficientes \(\beta_4\) y \(\beta_7\) tienen valores P bastante altos como para decir sin más que son significativos (13.16% y 11.04% respectivmente) pero sin embargo no son tan altos como para decir definitivamente que no son significativos. Cabe destacar que el \(\beta_4\) corresponde a la dummy D4 correspondida a la categoria Oeste de la variable region y a su vez el \(\beta_7\) corresponde a la variable (D4X1) que es la interacción con la variable ln(salario) que sucede cuando la categoria es Oeste. Por lo que, el hecho de que el cliente resida en la región Oeste (respecto de que resida en la categoría base Sur), podría ser significativo en la variable de respuesta ln(montoAnual), con 13.16% de riesgo de equivocarse al decir esto. Y tambien se puede decir que, que resida en la región Oeste, influirá en cuanto aumenta o disminuye la media del ln(montoAnual) por cada unidad que aumente o disminuya el ln(salario) de este cliente.* En un caso real, buscaría tomar más muestras de clientes de la región Oeste para poder aumentar la sensibilidad y la potencia del test de significación de \(\beta_4\) y \(\beta_7\), para poder llegar a una mejor conclusión.
En cuanto al resto de los coeficientes, todos tienen valores P demasiado altos como para decir que son significativos. Es decir, que es demasiado riesgoso rechazar la hipótesis nula. Por lo tanto, puedo decir que el resto de los coeficientes NO son significativos, es decir que las variables correspondidas NO aportan información suficiente al NO tener un grado de asociación lineal significativo con la variable ln(montoAnual).
Finalmente, por Principio de Parsimonia, dejo en el modelo aquellas variables de las que puedo decir que presentan un grado de asociación lineal con ln(montoAnual) con una baja probabilidad de error de tipo I. El modelo queda:
\(Y_i=\beta_0+\beta_1*X_{1i}+\tilde{\epsilon_i}\)
La recta estimada queda:
\(\hat{Y_i}=1.2505+0.9222*X_{1i}\)
b) Variable Dependiente: ln(montoAnual) Independientes: ln(salario) y frecuencia.
library(readxl)
Parte_2_B <- read_excel("C:/Users/Ignac/Desktop/Estadistica Aplicada II/TPs/TP2/Parte II b (usar este para R).xlsx")
head(Parte_2_B)Las variables ln(montoAnual) y ln(salario) van a estar definidas de igual manera que en el inciso anterior.
Defino la variable restante:
\(X_{2i}\):(frecuencia) frecuencia de compra del i-ésimo cliente [CUALITATIVA]. Variable considerada determinística, explicativa.
Como en el inciso anterior, aplico logaritmo a las observaciones de montoAnual y de salario para convertirlas a ln(montoAnual) y ln(salario):
library(dplyr)
Parte_2_B_ln <- transmute(Parte_2_B, ln_salario = log(salario), frecuencia = frecuencia, ln_montoAnual = log(montoAnual))
head(Parte_2_B_ln) # visualizo las primeras 6 filas de los datos habiendo aplicado el logaritmo naturalArmo un gráfico de dispersión de los datos, separando por colores y formas los puntos correspondientes a cada categoria de la variable frecuencia para un mejor entendimiento de los mismos:
library(ggplot2)
ggplot(Parte_2_B_ln, aes(x = ln_salario, y = ln_montoAnual, shape = frecuencia, color = frecuencia)) +
geom_point() +
labs(title="Diagrama de dispersión", subtitle="con la variable categórica Frecuencia",
x="ln(Salario del cliente)", y="ln(Monto anual gastado por el cliente)", shape="Frecuencia", colour="Frecuencia")
A simple vista, parece haber una tendencia a que los clientes con más alto nivel de salario y más monto anual gastado, compran con alta frecuencia. A su vez se observa que sucede algo parecido para los clientes que compran habitualmente (niveles medios de monto gastado y salario) y esporadicamente (niveles bajos). Aunque no es tan marcado como el caso anterior, ya que por ejemplo la nube de triangulos verdes (Esporádico) se llega a extender a niveles “altos” de monto anual gastado y de salario de clientes.
Codificación de la variable categórica, especificación del modelo de RLM e interpretación de los coeficientes de regresión:
Teniendo una variable cualitativa, frecuencia, la cual define un ordenamiento y por ende se la puede considerar como una variable ordinal. Con esto voy a considerar que la “distancia” entre cada categoría de la variable es equidistante. Entonces, codifico la variable cualitativa como una cuantitativa, representando los tres niveles del factor cualitativo frecuencia con los siguientes factores cuantitativos, definiendo la variable ordinal \(X_2\):
ordinal_B <- matrix(data = c(1,2,3), ncol=1, nrow=3,
dimnames=list(c("Esporádico", "Habitual", "Alta Frecuencia"), c("X2")))
print(ordinal_B)## X2
## Esporádico 1
## Habitual 2
## Alta Frecuencia 3Entonces planteo la especificación del modelo de RLM con esta variable ordinal del siguiente modo:
\(Y_i=\beta_0+\beta_1*X_{1i}+\beta_2*X_{2i}+\tilde{\epsilon_i}\)
Esto cuantifica el efecto sobre la ordenada al origen (efecto aditivo) que genera la variable categórica frecuencia modelada como una variable ordinal. La interpretación del coeficiente \(\beta_2\) es muy parecida a la interpretación de los coeficientes del inciso anterior \(\beta_k\). Es decir, este coeficiente se puede interpretar como que, independientemente del valor de la variable ln(salario) (es decir, independientemente del salario del cliente) (\(X_{1i}\)), la media de la variable ln(montoAnual) \(Y_i\) se va a incrementar en \(\beta_2*X_2\) con X2 tomando el valor respectivo a su categoría, es decir que la ordenada al origen quedaría \((\beta_0+\beta_2*X_2)\).
Para que se entienda mejor esto, expreso como quedaría la media de la variable ln(montoAnual) para cada categoría posible de la variable frecuencia:
Si la categoria es Esporádico: \(E[Y|X_2=1]=(\beta_0+\beta_2)+\beta_1*X_1\)
Si la categoria es Habitual: \(E[Y|X_2=2]=(\beta_0+2*\beta_2)+\beta_1*X_1\)
Si la categoria es Alta Frecuencia: \(E[Y|X_2=3]=(\beta_0+3*\beta_2)+\beta_1*X_1\)
En cuanto a la interpretación del coeficiente \(\beta_1\), correspondiente a ln(salario), esta va a ser la misma que ya se explicó en el inciso anterior.
Con la variable ordinal el dataframe queda así:
frecuencia <- Parte_2_B_ln$frecuencia
X2 <- c()
for (i in 1:length(frecuencia)) {
if (frecuencia[i] == "Esporadico") {
X2 <- c(X2, 1)
}
if (frecuencia[i] == "Habitual") {
X2 <- c(X2, 2)
}
if (frecuencia[i] == "Alta Frecuencia") {
X2 <- c(X2, 3)
}
}
Parte_2_B_ln_ordinal <- data.frame("ln_salario"=Parte_2_B_ln$ln_salario,
"frecuencia"=frecuencia,"ln_montoAnual"=Parte_2_B_ln$ln_montoAnual,
"X2"=X2)
head(Parte_2_B_ln_ordinal)Saco la variable categórica frecuencia porque ya no es necesaria:
Parte_2_B_ln_ordinal <- data.frame("ln_salario"=ln_salario, "ln_montoAnual"=ln_montoAnual,
"X2"=X2)Ahora calculo los coeficientes de regresión con los cálculos matriciales:
# Armo el VECTOR DE OBSERVACIONES DE Y
ln_montoAnual_B <- Parte_2_B_ln_ordinal$ln_montoAnual
ln_montoAnual_B <- matrix(data = ln_montoAnual, nrow = length(ln_montoAnual))
# Armo la MATRIZ DE DATOS de dimensión (n x p)
n <- 100 # cantidad de obs de cada variables
p <- 3 # cantidad de parámetros a estimar
unos_B <- c()
for (i in 1:n) {
unos_B <- c(unos_B, 1) # creo la columna de unos
}
ln_salario_B <- Parte_2_B_ln_ordinal$ln_salario # columna de ln(salario)
matriz_datos_B <- matrix(data = c(unos_B, ln_salario_B, X2), nrow=n, ncol=p)
head(matriz_datos_B) # muestro las primeras 6 filas de la matriz de datos## [,1] [,2] [,3]
## [1,] 1 9.532424 1
## [2,] 1 11.247657 2
## [3,] 1 10.515967 1
## [4,] 1 11.415313 3
## [5,] 1 11.513925 3
## [6,] 1 10.599132 2Ahora armo la matriz \((X^{T}*X)^{-1}\), nuevamente la “matriz normal”:
matriz_normal_B <- solve(t(matriz_datos_B) %*% matriz_datos_B)Calculo el Vector de Coeficientes Estimado \(\beta=(X^{T}*X)^{-1}*X^T*Y\).
vect_coef_estimado_B <- matriz_normal_B %*% t(matriz_datos_B) %*% ln_montoAnual_B
rownames(vect_coef_estimado_B) <- c("Beta0", "Beta1", "Beta2")
colnames(vect_coef_estimado_B) <- c("Estimación Puntual")
vect_coef_estimado_B## Estimación Puntual
## Beta0 5.4400709
## Beta1 0.4140851
## Beta2 0.7052814Tests de significación para los coeficientes y elección del modelo
Se realizará inferencia sobre los coeficientes estimados, usando el mismo estadístico que se usó en el inciso anterior (t de Student).
Test de hipótesis UNIlateral para el coeficiente \(\beta_1\):
\(H_0)\beta_1 \leq 0\) ; \(H_a)\beta_1 > 0\)
Dado que se puede suponer extra-estadísticamente que el ln(montoAnual) debería aumentar si aumenta ln(salario) entonces la hipótesis alternativa para este coeficiente será que este tenga que ser mayor que cero (pendiente positiva).
Test de hipótesis BIlateral para los otros 2 coeficientes estimados:
\(H_0)\beta_j=0\) ; \(H_a)\beta_j \neq 0\)
\((j \neq 1)\)
Calculo y armo una tabla similar a la del inciso anterior para calcular los valores P:
predichos_B <- matriz_datos_B %*% vect_coef_estimado_B # predichos para los datos observados
Q_B <- sum((ln_montoAnual_B-predichos_B)^{2}) # suma de cuadrado de los residuos
S_B <- sqrt(Q_B/(n-p)) # desvío residual
# Calculo el desvío de cada coeficiente estimado
C_ii_B <- diag(matriz_normal_B) # elementos de la diagonal de la "matriz normal"
desvio_coef_estimado_B <- S_B*sqrt(C_ii_B)
# Calculo los valores p (los calculo todos igual y despues modifico el de beta1)
vect_tobs_B <- (vect_coef_estimado_B-0)/desvio_coef_estimado_B # estadísticos t observados
valoresP_B <- 2*pt(abs(vect_tobs_B), df=n-p, lower.tail=FALSE) # valores P test BIlateral
valoresP_B[2] = valoresP_B[2]/2 # modifico el valor p para el beta 1 (test unilateral)
# Calculo los IC del 95% (NC=0.95=1-alfa=1-0.05)
LI_coef_estimado_B <- vect_coef_estimado_B-qt(1-0.05/2, df=n-p)*desvio_coef_estimado_B
LS_coef_estimado_B <- vect_coef_estimado_B+qt(1-0.05/2, df=n-p)*desvio_coef_estimado_B
# resultados de los test de significación
resultados_B <- c("---", "SIGNIFICATIVO", "SIGNIFICATIVO")
matrix(data = c(vect_coef_estimado_B, desvio_coef_estimado_B, LI_coef_estimado_B, LS_coef_estimado_B, vect_tobs_B, valoresP_B, resultados_B),
nrow = 3, ncol = 7,
dimnames = list(c("Beta0", "Beta1", "Beta2"),
c("Estimación puntual", "Desvío estimado del coeficiente", "LI (IC95%)", "LS (IC95%)", "t observado", "Valor P", "Resultado")))## Estimación puntual Desvío estimado del coeficiente LI (IC95%)
## Beta0 "5.44007091602941" "1.09256156138568" "3.27163865299961"
## Beta1 "0.414085126546392" "0.112709336886565" "0.190388292347365"
## Beta2 "0.705281381330588" "0.0877461949578861" "0.531129473713168"
## LS (IC95%) t observado Valor P
## Beta0 "7.60850317905921" "4.97918937321011" "0.00000277455973159209"
## Beta1 "0.637781960745419" "3.67392035109871" "0.000195654011620172"
## Beta2 "0.879433288948008" "8.0377431940962" "0.00000000000223100487469863"
## Resultado
## Beta0 "---"
## Beta1 "SIGNIFICATIVO"
## Beta2 "SIGNIFICATIVO"NOTA: el valor P correspondiente al \(\beta_1\) es para test UNIlateral, fue modificado dado \(\alpha^{\star}_{unilat}=\frac{\alpha^{\star}_{bilat}}{2}\). El resto de los valores P sin para tests BIlaterales.
Para este caso todos los coeficientes dan significativos ya que tienen valores P muy bajo, aunque cabe destacar que el valor P para el test del coeficiente \(\beta_2\) es mucho menor que el del \(\beta_1\). Por lo tanto, puedo decir que tanto la variable ln(salario) \(X_1\), como la variable ordinal \(X_2\) (que representa a la categórica frecuencia), tienen un cierto grado de asociación lineal con la variable de respuesta ln(montoAnual), con solo un 0.01957% y 0.0000000002231% de probabilidad de cometer error de tipo I al decir esto, respectivamente.
Por ende, para el modelo final no voy a sacar ninguna de las variables porque considero que todas son significativas y aportan información útil para la regresión:
\(Y_i=\beta_0+\beta_1*X_{1i}+\beta_2*X_{2i}+\tilde{\epsilon_i}\)
La recta estimada queda:
\(\hat{Y_i}=5.4401+0.4141*X_{1i}+0.7053*X_{2i}\)
Parte III
Se debe utilizar el modelo completo, es decir incluyendo todas las variables asignadas en la planilla Excel para la Parte III. No debe quitar ni agregar ninguna variable.
library(readxl)
Parte_3 <- read_excel("C:/Users/Ignac/Desktop/Estadistica Aplicada II/TPs/TP2/Parte III.xlsx")
head(Parte_3)Defino las variables del modelo completo:
\(Y_i\): Monto gastado en compras, por el i-ésimo cliente en el último año (montoAnual)[$]. Variable aleatoria, de respuesta/explicada.
\(X_{1i}\): Salario del i-ésimo cliente (salario)[$]. Variable considerada determinística, explicativa.
\(X_{2i}\): Edad del i-ésimo cliente (edad)[años]. Variable considerada determinística, explicativa.
\(X_{3i}\): Cantidad de catálogos promocionales enviados durante los últimos 2 años al i-ésimo cliente (nCatalogos). Variable considerada determinística, explicativa.
\(X_{4i}\): Indica si el i-ésimo cliente vive en el barrio (SI/NO) (cerca)[CUALITATIVA]. Variable considerada determinística, explicativa.
La especificación de este modelo de RLM, llamado “completo”, es la siguiente:
\(Y_i=\beta_0+\beta_1*X_{1i}+\beta_2*X_{2i}+\beta_3*X_{3i}+\beta_4*X_{4i}+\tilde{\epsilon_i}\)
a) Estime el modelo completo
Teniendo una variable cualitativa, cerca, con dos categorias posibles (Si o No), armo una dummy con categoría base No para codificarla, la llamo D4:
cerca <- Parte_3$cerca
de_4 <- c()
for (i in 1:length(cerca)) {
ifelse(cerca[i] == "Si", # la D4 corresponde a la categoría Si
de_4 <- c(de_4, 1), de_4 <- c(de_4, 0))
}
Parte_3 <- data.frame("montoAnual"=Parte_3$montoAnual, "salario"=Parte_3$salario, "edad"=Parte_3$edad,
"nCatalogos"=Parte_3$nCatalogos, "cerca"=de_4)
head(Parte_3)Especificación del modelo, recta de regresión poblacional:
\(Y_i=\beta_0+\beta_1*X_{1i}+\beta_2*X_{2i}+\beta_3*X_{3i}+\beta_4*D_{4i}+\tilde{\epsilon_i}\)
La recta de regresión estimada es entonces:
\(\hat{Y_i}=b_0+b_1*X_{1i}+b_2*X_{2i}+b_3*X_{3i}+b_4*D_{4i}\)
Calculo los coeficientes de regresión estimados matricialmente con el método de mínimos cuadrados:
options(scipen = 999)
# VECTOR DE OBSERVACIONES Y
observaciones <- Parte_3$montoAnual
observaciones <- matrix(data = observaciones, nrow = length(observaciones), ncol = 1)
# MATRIZ DE DATOS
filas_n <- 100 # cantidad de obs de cada variables
columnas_p <- 5 # cantidad de parámetros a estimar
vector_unos <- c()
for (i in 1:filas_n) {
vector_unos <- c(vector_unos, 1) # creo la columna de unos
}
equis_1 <- Parte_3$salario
equis_2 <- Parte_3$edad
equis_3 <- Parte_3$nCatalogos
m_datos <- matrix(data = c(vector_unos, equis_1, equis_2, equis_3, de_4),
nrow=filas_n, ncol=columnas_p) # matriz de datos
# MATRIZ (Xt*X)^{-1} (la llamo "matriz normal")
m_normal <- solve(t(m_datos) %*% m_datos)
# COEFICIENTES ESTIMADOS
estimaciones <- m_normal %*% t(m_datos) %*% observaciones
rownames(estimaciones) <- c("Beta0", "Beta1", "Beta2", "Beta3", "Beta4")
colnames(estimaciones) <- c("Estimación Puntual")
estimaciones## Estimación Puntual
## Beta0 -108538.049115
## Beta1 2.500079
## Beta2 1129.308599
## Beta3 5913.708705
## Beta4 -77517.149565La recta estimada queda:
\(\hat{Y_i}=-108538.05+2.5*X_{1i}+1129.31*X_{2i}+5913.71*X_{3i}-77517.15*D_{4i}\)
b) Analice la colinealidad del modelo. ¿Haría alguna recomendación al respecto?
Análisis de Multicolinealidad:
Deseo analizar si hay prescencia de asociación lineal entre las variables explicativas del modelo, para eso calculo la matriz de correlaciones de las variables explicativas (para esto utilizo funciones de R):
explicativas_parte_3 <- data.frame("salario"=Parte_3$salario, "edad"=Parte_3$edad,
"nCatalogos"=Parte_3$nCatalogos, "cerca"=de_4)
m_corr <- round(cor(x = explicativas_parte_3, method = "pearson"), 3)
m_corr## salario edad nCatalogos cerca
## salario 1.000 0.111 0.235 -0.058
## edad 0.111 1.000 -0.043 -0.001
## nCatalogos 0.235 -0.043 1.000 -0.094
## cerca -0.058 -0.001 -0.094 1.000Esta matriz me da los coeficientes de correlación entre cada variable en las filas y cada variable en las columnas, por eso es simétrica respecto de la diagonal.
Ahora armo un gráfico que me permite ver tanto estos coeficientes de correlación entre las variables como así tambien los diagramas de dispersión entre las variables, lo cual permite una ver gráficamente (ademas de cuantitativamente) si hay asocianción lineal entre variables explicativas. Ademas, este gráfico da los histogramas de cada variable.
library(ggplot2)
library(GGally)## Warning: package 'GGally' was built under R version 4.0.3## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2ggpairs(explicativas_parte_3,
diag = list(continuous = "barDiag"), axisLabels = "none")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Los valores de los coeficientes de en este gráfico son bastante bajos, y tampoco se observan en los gráficos que parezca haber una asociación lineal entre estas variables explicativas. El mayor coeficiente de correlación lineal que se observa es entre las variables salario y nCatalogos, que es de 0.235 (marcado con un asterísco en el gráfico), pero sin embargo este no pareciera ser una correlación muy alta.
Para calcular un buen indicador global de multicolinealidad a partir de esta matriz, calculo el DET. Es es el determinante de la matriz de correlaciones.
det(m_corr)## [1] 0.9188914\(DET=0.9188914\)
Este determinante es muy cercano a 1, por lo que es muy baja la colinealidad entre estas variables, globalmente hablando. Se considera que si \(DET \geq 0.2\) entonces NO hay multicolinealidad, y este DET es por mucho mayor a 0.2.
Ahora, aunque ya se ve que no hay una gran multicolinealidad, puedo ver mejor cuales son las variables que traen esta multicolinealidad con el cálculo de los VIF (factores de inflación de varianza):
\(VIF_j=\frac{1}{1-R^{2}_j}\)
Siendo \(R^{2}_j\) el coeficiente de determinación donde Xj es la variable de respuesta y el resto las explicativas.
Este VIF me indica cuánto se multiplica la varianza del coeficiente \(beta_j\), en comparación a si \(\beta_j=0\).
Calculo los VIF, para mayor facilidad uso las funciones de R, principalmente voy a usar la funcion de modelo lineal (lm) que da toda la información necesaria (coeficientes, predichos, residuos, etc):
modelo <- lm(formula = montoAnual ~ equis_1 + equis_2 + equis_3 + de_4, data = Parte_3)
library(car)## Warning: package 'car' was built under R version 4.0.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.0.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodevif(modelo)## equis_1 equis_2 equis_3 de_4
## 1.076731 1.017589 1.071542 1.010326Siendo que todos los VIF (por mucho) menores a 5, se dice que la multicolinealidad no es lo suficientemente elevada para afectar las estimaciones. Tal cual se habia previsto anteriormente. Se podría decir que las variables explicativas que traen la (muy pequeña) multicolinealidad son la X1 (salario) y la X3 (nCatalogos) porque tienen los VIF mas altos.
Por lo tanto, NO haria ninguna recomendación al respecto ya que la multicolinealidad entre las explicativas es muy baja.
Si hubiera multicolinealidad habría que analizar si este problema es muestral, cambiando la muestra o agregando observaciones que no esten muy correlacionadas con las anteriores, o si es poblacional, es decir que la dependencia es entre las explicativas es teórica y mi modelo funcionaria para esta población (no podria extrapolar a otras poblaciones parecidas). En este último caso debería observar los leverages detenidamente.
c) Realice un diagnóstico completo del modelo analizando:
a. Comportamiento de los residuos
Análisis de Residuos
Calculo y grafico los residuos estudentizados:
\(r_i=\frac{Y_i-\hat{Y_i}}{S*\sqrt{1-h_{ii}}}\)
En el numerador tengo los residuos ei y en el denominador esta la estimación del desvío de los residuos, es como si estuviera “estandarizando” con la estimación del desvio de ei por eso digo que son estudentizados. La S es el desvío residual y la hii son los leverage que es la diagonal de la matriz H que calculo explicitamente mas adelante cuando analice los leverages.
Estos residuos son mejores que los ordinarios o los estandarizados. Ya que permiten un mejor análisis porque aumentan su valor cuanto más alejados estén del centroide de las variables explicativas (mayor influencia).
Calculo y grafico los residuos estudentizados estudentizados vs los predichos:
r_estudentizados <- rstudent(modelo) # aca calculo los residuos estudentizados
ggplot(data = Parte_3, aes(modelo$fitted.values, r_estudentizados)) +
geom_point() +
geom_smooth(color = "blue", se = FALSE) + # promedio de los residuos
geom_hline(yintercept = 0) +
geom_hline(yintercept = 2, color = "red") +
geom_hline(yintercept = -2, color = "red") +
labs(title="Residuos Estudentizados VS Predichos", x="Predichos", y="Residuos Estudentizados")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
La línea azúl indica (que esta estimada por mínimos cuadrados, pero no tiene que ver con una regresión) sirve para observar que los residuos estudentizados presentan una forma bastante curvilínea, con lo cual debería probarse aplicar una transformación lineal para corregir este incombeniente. Esto indica que claramente no están variando aleatoriamente respecto del valor nulo, como deberían según el supuesto de Auscencia de Sesgo, y además se ve que aumenta mucho la dispersión de los puntos a mayores valores de los predichos, por lo que tampoco se cumple el supuesto de Homocedasticidad de los residuos.
b. Análisis de Outliers
Lo primero para el análisis de outliers es observar el gráfico que hice anteriormente, en el que puse dos líneas rojas (en los valores de residuos estudentizados 2 y -2) que indican los límites que se consideran “normales” para los residuos estudentizados. Por fuera de estos límites los puntos son posibles outliers, es decir, que corresponden a valores de la variable de respuesta Y que son condicionalmente inusuales dados los valores de las variables explicativas. El problema con estos es que pueden afectar seriamente a las estimaciones de la variable de respuesta al afectar a las estimaciones de los parámetros, porque podría cambiar la inclinación del hiperplano estimado de regresión múltiple. Estos puntos corresponden a las observaciones 47, 30, 8, 88 y 84 (habia puesto los números de observaciones en el gráfico para fijarme esto pero los saqué porque sino no se aprecia bien). Se puede decir que presentan una elevada discrepancia. Esto no quiere decir que sean definitivamente outliers, sino que son candidatos a serlo. En los proximos incisos voy a definir si lo son o no, cuantitativamente.
c. Análisis de leverage e influencia
Primero hago un análisis de leverage o balanceo, estos son la distancia que hay desde el centroide de las variables explicativas hasta el punto (observación) que se esté analizando. Para RLM, los leverages van a ser la diagonal (hii) de la matriz H, la cual proyecta los valores observados de Y en valores pronosticados y se calcula solo operando la matriz de datos.
\(H=X*(X^{T}*X)^{-1}*X^{T}\)
=> \(\hat{Y}=H*Y\)
El valor del leverage va a ser un importante indicador de si el punto es influyente y por ende tambien si es un outlier o no. Para saber si el leverage de un punto es alto o no, se usa el siguiente criterio:
El punto tiene un alto leverage si \(h_{ii} \geq 2*\frac{p}{n}\)
Entonces, armo un gráfico de los leverages (vs el número de observación) y marco con una línea roja el límite que da el criterio:
leverage <- hatvalues(modelo)
caso <- seq(1, 100, 1)
leve_casos_parte_3 <- data.frame("caso"=caso, "leverage"=leverage)
ggplot(leve_casos_parte_3, aes(x=caso, y=leverage)) +
geom_point() +
geom_text(
label=rownames(leve_casos_parte_3),
check_overlap = TRUE) +
geom_hline(yintercept=2*(columnas_p/filas_n), color="red") + #VALOR LÍMITE
labs(title = "Leverage o Balanceo")
El gráfico muestra que no muestra que haya puntos con suficiente leverage como para decir que su leverage es alto según el criterio. Pero sin embargo, se observa que los de las observaciones/casos 37 y 84 tienen altos leverages relativamente con el resto y son muy cercanos al límite de ser altos absolutamente. Por lo que serán puntos a tener en cuenta en el siguiente análisis.
Para poder hacer un análisis de influencia, hay que tener en cuenta que, el hecho de que un punto sea influyente quiere decir que afecta severamente a la estimación de los coeficientes. Y por ende para definir la influencia de un punto hay que tener en cuenta tanto su discrepancia (valor de los residuos) como su leverage.
influencia = leverage * discrepancia
Para medir cuantitativamente la influencia de un punto se puede utilizar la distancia de Cook, la cual mide la influencia en los valores predichos sise elimina la observación analizada. Es decir que va a ser un indicador global de influencia, teniendo en cuenta tanto el tamaño del residuo (discrepancia) como el leverage.
\(D_i=\frac{(\hat{Y_i}-\hat{Y_{-i})^2}}{p*S^{2}}=\frac{r_i^{2}}{k+1}*\frac{h_{ii}}{1-h_{ii}}\)
En esta expresión se ve lo dicho anteriormente de que influencia = leverage * discrepancia.
Por lo tanto calculo y grafíco las distancias de Cook para cada punto:
D_cook <- cooks.distance(modelo)
cook_casos_parte_3 <- data.frame("caso"=caso, "distancia de Cook"=D_cook)
ggplot(cook_casos_parte_3, aes(x=caso, y=D_cook)) +
geom_point() +
geom_text(
label=rownames(cook_casos_parte_3),
check_overlap = TRUE) +
labs(title = "Distancia de Cook", y="Distancia de Cook")
Primero tengo en cuenta que, según el criterio general, se deben analizar las distancias de cook mayores a la unidad, y tal cuan se ve en este diagrama, todas estan muy por debajo de la unidad. Sin embargo, tambien hay que observar si hay algún punto que tenga una distancia de Cook mucho mayor al resto, porque quiere decir que ese punto será mucho más influyente, y esto se puede ver en que sucede para el caso 84. Y hay que destacar que este punto 84 dio que tiene tanto alto leverage como un alto valor de su residuo correspondiente. Entonces debo hacer un analisis más profundo para este punto.
Para ver si verdaderamente la observación 84 es influyente, voy a los DFBETAS y los DFFITS para este punto.
Los DFBETAS miden la influencia en los coeficientes estimados si se elimina la observación analizada:
\(DFBETAS_{ij}=\frac{b_j-b_{j(-i)}}{S_{b_{j(-i)}}}\)
Criterio de punto potencialmente influyente:
\(|DFBETAS_{ij}|>\frac{2}{\sqrt{n}}\)
print(c("Límite:", 2/sqrt(filas_n)))## [1] "Límite:" "0.2"Con las funciones de R puedo calcular los DFBETAS con facilidad, pero en este caso solo analizo la observación 84:
DFBETAS <- dfbetas(modelo)
DFBETAS[84,]## (Intercept) equis_1 equis_2 equis_3 de_4
## -0.6683906 0.3084594 0.6955104 0.4153536 -0.3683626Los DFFITS miden la influencia en los valores predichos si se elimina la observación analizada:
\(DFFITS_{i}=\frac{\hat{Y_{i}}-\hat{Y_{-i}}}{S_{-i}*\sqrt{h_{ii}}}\)
Criterio de punto potencialmente influyente:
\(|DFFITS_{i}|>2*\sqrt{\frac{p}{n}}\)
print(c("Límite:", round(2*sqrt(columnas_p/filas_n), 3)))## [1] "Límite:" "0.447"Con las funciones de R puedo calcular los DFBETAS con facilidad, pero en este caso solo analizo la observación 84:
DFFITS <- dffits(modelo)
DFFITS[84]## 84
## 1.080646Finalmente se tiene que, para la observación 84, los DFBETAS para todos los coeficientes tienen módulos mayores al límite 0.2. Y que el DFFITS es tambien mayor al límite.
Por lo tanto con todas estas mediciones es muy probable que la observación 84 sea un punto influyente.
d. Cumplimiento de los supuestos
Supuesto de Ausencia de Vicio (insesgadez) Como ya se habló, se ve en el diagrama de residuos estudentizados vs predichos que este supuesto no se cumple. Debido a que el gráfico muestra que los puntos forman un patron geométrico (curva), es decir, que hay un sesgo sistemático de los residuos. Y como estos son estimaciones de las perturbaciones, entonces no se cumple que la media de las perturbaciones es nula.
Supuesto de Homocedasticidad (igualdad de varianza de perturbaciones) Como también se mencionó cuando se hizo el gráfico de residuos estudentizados, la dispersión de los puntos en el gráfico no es uniforme, ya que estos no varian aleatoriamente respecto del cero. Por lo tanto digo que este supuesto tampoco se cumple.
Supuesto de Normalidad de las Perturbaciones
Este supuesto plantea que \(\epsilon_i \sim Normal(0, \sigma^{2})\)
Para ver si se cumple voy a utilizar un gráfico qqplot, el cual compara valores teóricos de una distribución con los valores que se quiere saber si siguen esa distribución. En este caso quiero saber si los residuos siguen o no una distribución Normal, por ende los grafico del siguiente modo:
qqnorm(residuals(modelo), col="blue")
qqline(residuals(modelo), col="red")
Se ve que en general los residuos se ajustan bastante bien a la distribución Normal. Sin embargo, se observan en los valores mas altos de los residuos, que hay unos 4 puntos que se desvían de la distribución Normal teórica. Estos puntos de alto valor de residuo, muy probablemente sean los mismos posibles outliers de los que se habló en el gráfico de resifuos estudentizados vs predichos. Entonces, no se puede decir que este supuesto no se cumple, aunque sin embargo estos puntos pueden generar un ruido inesperado en las estimaciones, lo cual es problemático.
Supuesto de Ausencia de Autocorrelación
Este supuesto implica que \(Cov(\epsilon_i;\epsilon_j)=0\), es decir, que no hay autocorrelación entre las perturbaciones. Entonces tampoco deberia haberla entre los residuos. Suele suceder que sí hay autocorrelación si los datos están ordenados en el tiempo o espacio. Para este caso de aplicación no parecería que suceda esto. Si se quiere saber cuantitativamente si se cumple o no el supuesto, se podría hacer un test para el coeficiente de correlación, que en la hipótesis nula diga que es cero. Y se usaría el estadístico de Durbin-Watson, pero en esta materia no vemos como hacer esto.
d) Para la última observación del set de datos:
a. Calcule el DIFFITS y los DIFBETAS con Excel. ¿Estos valores le parecen elevados o no? ¿Por qué? Interprete estos valores en términos del problema.
ANEXO DE EXCEL CON LOS CÁLCULOS
De los calculos de excel para la observación número 100, tengo que:
\(DFBETAS_{100;beta0}=-0.01944\)
\(DFBETAS_{100;beta1}=-0.03671\)
\(DFBETAS_{100;beta2}=0.07236\)
\(DFBETAS_{100;beta3}=-0.03013\)
\(DFBETAS_{100;beta4}=0.02142\)
\(DFFITS_{100}=0.09999\)
como ya se mencionó, los criterios para saber si los valores son elevados o no son los siguientes:
\(|DFBETAS_{ij}|>(\frac{2}{\sqrt{n}}=0.2)\)
\(|DFFITS_{i}|>(2*\sqrt{\frac{p}{n}}=0.447)\)
Por lo tanto ni los DFBETAS ni el DFFITS presentan valores elevados. Esto es porque parece ser que la influencia de la observación 100 es baja. Ya que, para esta observación, los DFBETAS son una medición de la influencia en los coeficientes estimados si se elimina esta observación 100, osea que calculan cuanto cambian las pendientes del hiperplano o la ordenada al origen, si no considero el punto. Y el DFFITS es la medición de la influencia en los valores predichos si se elimina la observación 100, es decir, es un valor de cuánto cambia el pronóstico, medido en cantidades de desvios estándar, si no considero el punto. Entonces, sin la observación 100, el pronóstico me va a cambiar en 0.447 desvíos estándar.
e) Para el individuo cuyos datos se reproducen en la hoja “ParteIIIpron”, indique entre valores se encontrarán sus compras a la empresa con 90% de confianza. ¿Considera que se extrapolando al realizar esta estimación?
Datos a pronosticar:
X_0 <- c("salario"=140700, "edad"=32, "nCatalogos"=12, "cerca"=1)
X_0## salario edad nCatalogos cerca
## 140700 32 12 1X_0 <- matrix(data = c(1, X_0), nrow = 5)Armo un intervalo de predicción del 90%, NC = 0.9 = 1-0.1.
Inferencia sobre: \(Y|X_0 =\beta_0+\beta_1*X_{1;0}+\beta_2*X_{2;0}+\beta_3*X_{3;0}+\beta_4*X_{4;0}\)
Estadístico:
\(t_(n-p)=\frac{\hat{Y}|X_0-Y|X_0}{\sqrt{\hat{D}^{2}(\hat{Y_0})+S^2}}\)
Estimación puntual del montoAnual con los datos dados:
Ypico_dado_X0 <- t(X_0) %*% estimaciones
Ypico_dado_X0 # estimación puntual del montoAnual dados los datos## Estimación Puntual
## [1,] 272808.4Varianza residual S^2:
SC_residual <- sum((residuals(modelo))^{2})
varianza_residual <- SC_residual/(filas_n-columnas_p) #varianza residual
varianza_residual## [1] 5677309555varianza estimada de la estimación puntual del montoAnual:
Dpico_cuadrado_Ypico <- (t(X_0) %*% m_normal %*% X_0)*varianza_residual
Dpico_cuadrado_Ypico## [,1]
## [1,] 554890636Armo el intervalo de predicción del 90%:
LI_prediccion <-Ypico_dado_X0-qt(1-0.1/2, df=filas_n-columnas_p)*sqrt(Dpico_cuadrado_Ypico+varianza_residual)
LS_prediccion <-Ypico_dado_X0+qt(1-0.1/2, df=filas_n-columnas_p)*sqrt(Dpico_cuadrado_Ypico+varianza_residual)
c("LI_90%"=LI_prediccion , "LS_90%"=LS_prediccion )## LI_90% LS_90%
## 141677.8 403938.9Finalmente, el intervalo de predicción del montoAnual, para los datos dados es:
\(IP(Y|X_0)_{90\%}=[141677.8 , 403938.9]\)
Considero que sí se está extrapolando al realizar esta estimación, ya que los datos de la muestra, con la que se armó la estimación de los coeficientes, tiene valores de la variable salario que NO pasan los 140000. Y por ende como en los datos dados la variable salario es de 140700 sí se está extrapolando, aunque no por mucho, lo cual puede ser aceptable en determinados casos. Sin embargo, como extrapolo, me alejo del centroide de las explicativas y por ende se amplia mucho el intervalo de predicción (aún más que el de confianza). Por lo que esta estimación por intervalo de predicción es poco precisa, en comparación de si no extrapolaría.