1. Un economista del Departamento de Recursos Humanos de Florida State está preparando un estudio sobre el comportamiento del consumidor. Él recolectó los datos que aparecen en miles de dólares para determinar si existe una relación entre el ingreso del consumidor y los niveles de consumo. Los datos recolectados se pueden ver en la siguiente tabla:
ingreso <- c(24.3, 12.5, 31.2, 28.0, 35.1, 10.5, 23.2, 10.0, 8.5, 15.9, 14.7, 15.0)

consumo <- c(16.2, 8.5, 15.0, 17.0, 24.2, 11.2, 15.0, 7.1, 3.5, 11.5, 10.7, 9.2)

df <- data.frame(ingreso, consumo)
names(df) <- c("Ingreso anual","Consumo anual")
print(df)
##    Ingreso anual Consumo anual
## 1           24.3          16.2
## 2           12.5           8.5
## 3           31.2          15.0
## 4           28.0          17.0
## 5           35.1          24.2
## 6           10.5          11.2
## 7           23.2          15.0
## 8           10.0           7.1
## 9            8.5           3.5
## 10          15.9          11.5
## 11          14.7          10.7
## 12          15.0           9.2

El economísta se pregunta cuánto consumo anual tendrá una persona que tenga un ingreso anual de 29.7 mil dolares.

Primero, vamos a mostrar un gráfico de dispersión para ver si las variables pueden tener una distribución normal.

library(ggplot2)

ingreso <- c(24.3, 12.5, 31.2, 28.0, 35.1, 10.5, 23.2, 10.0, 8.5, 15.9, 14.7, 15.0)

consumo <- c(16.2, 8.5, 15.0, 17.0, 24.2, 11.2, 15.0, 7.1, 3.5, 11.5, 10.7, 9.2)

dfx <- data.frame(ingreso, consumo)

ggplot(data = dfx, aes(x = ingreso, y = consumo)) +
  geom_point(color = "blue", size = 3) +   # Puntos
  labs(title = "Ingreso vs Consumo", x = "Ingreso anual", y = "Consumo anual") +
  theme_minimal()   # Estilo del gráfico

Viendo el gráfico, podemos suponer que efectivamente, los datos tienen una distribución normal, sin embargo, ahora hacemos la prueba y de paso verémos si tienen correlación.

H_0: No existe una correlación entre el ingreso anual y el consumo anual.

H_A: Existe una correlación entre el ingreso anual y el consumo anual.

ingreso <- c(24.3, 12.5, 31.2, 28.0, 35.1, 10.5, 23.2, 10.0, 8.5, 15.9, 14.7, 15.0)

consumo <- c(16.2, 8.5, 15.0, 17.0, 24.2, 11.2, 15.0, 7.1, 3.5, 11.5, 10.7, 9.2)

shapiro.test(ingreso)
## 
##  Shapiro-Wilk normality test
## 
## data:  ingreso
## W = 0.91209, p-value = 0.2269
shapiro.test(consumo)
## 
##  Shapiro-Wilk normality test
## 
## data:  consumo
## W = 0.96663, p-value = 0.8725
cor.test(ingreso, consumo, method = "pearson") #Ya que ambas variables se distribuyen normalmente
## 
##  Pearson's product-moment correlation
## 
## data:  ingreso and consumo
## t = 7.3762, df = 10, p-value = 2.38e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7305293 0.9774317
## sample estimates:
##       cor 
## 0.9190975

Como el p-valor es menor que 0.05, entonces, rechazamos nuestra H_0 y aceptamos H_A, por lo tanto, existe una correlación entre el ingreso anual y el consumo anual. Además, como cor > 0.8 decimos que tiene una correlación positiva fuerte, esto significa que, a mayor ingreso anual, mayor gasto anual.

Ahora, harémos un modelo de regresón lineal para intentar predecir el gasto anual de una persona sabiendo su ingreso anual.

ingreso <- c(24.3, 12.5, 31.2, 28.0, 35.1, 10.5, 23.2, 10.0, 8.5, 15.9, 14.7, 15.0)

consumo <- c(16.2, 8.5, 15.0, 17.0, 24.2, 11.2, 15.0, 7.1, 3.5, 11.5, 10.7, 9.2)

dfx <- data.frame(ingreso, consumo)

modelo <- lm(consumo ~ ingreso, data = dfx)

summary(modelo)
## 
## Call:
## lm(formula = consumo ~ ingreso, data = dfx)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1928 -0.5426  0.0088  0.8500  3.5613 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.77788    1.58292   1.123    0.288    
## ingreso      0.55817    0.07567   7.376 2.38e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.251 on 10 degrees of freedom
## Multiple R-squared:  0.8447, Adjusted R-squared:  0.8292 
## F-statistic: 54.41 on 1 and 10 DF,  p-value: 2.38e-05

Donde el Estimate es el valor esperado del consumo anual cuando el ingreso anual es 0 y por cada aumento de 1 unidad en el ingreso anual, el consumo anual aumenta en 0.55817, en otras palabras, nuestro modelo es: f(x) = 0.55817x + 1.77788. Y la fiabilidad de este modelo es del 84.47% (R-squared). Vamos a verlo de una manera más gráfica:

library(ggplot2)

ingreso <- c(24.3, 12.5, 31.2, 28.0, 35.1, 10.5, 23.2, 10.0, 8.5, 15.9, 14.7, 15.0)

consumo <- c(16.2, 8.5, 15.0, 17.0, 24.2, 11.2, 15.0, 7.1, 3.5, 11.5, 10.7, 9.2)

dfx <- data.frame(ingreso, consumo)


ggplot(dfx, aes(x = ingreso, y = consumo)) +
  geom_point(color = "blue", size = 3) +           
  geom_abline(
    intercept = 1.77788, slope = 0.55817,            # Función f(x) = 0.542x + 1.633
    color = "red", linewidth = 1, linetype = "solid"
  ) +
  labs(
    title = "Regresión lineal: Consumo ~ Ingreso",
    x = "Ingreso (miles de dólares)",
    y = "Consumo (miles de dólares)"
  ) +
  theme_minimal()

Por último, para responder la preguntada planteada, simplemente reemplazamos en x = 29.7

x = 29.7

 0.55817*x + 1.77788
## [1] 18.35553

Por lo tanto, una persona que tenga un ingreso de 29.7 mil dolares anuales, tendrá un gasto aproximado de 18.36 mil dolares anuales.

  1. Regresión Theil-Sen.
set.seed(123) 


vector1 <- rexp(100, rate = 0.5)


vector2 <- 2.5 * vector1 + rchisq(100, df = 3)


df2 <- data.frame(vector1, vector2)

print(df2)
##          vector1   vector2
## 1    1.686914522  9.012492
## 2    1.153220542  5.848106
## 3    2.658109736  8.046527
## 4    0.063154718  2.591927
## 5    0.112421952  1.422732
## 6    0.633002433  4.088466
## 7    0.628454584  2.512457
## 8    0.290533608  6.145973
## 9    5.452472929 14.811512
## 10   0.058306894  6.385436
## 11   2.009660115  6.580527
## 12   0.960429455  4.785299
## 13   0.562027255  8.544512
## 14   0.754235662  3.961665
## 15   0.376568082  5.985607
## 16   1.699572259  5.287840
## 17   3.126407079  9.726976
## 18   0.957520833  5.400028
## 19   1.181869671  6.147174
## 20   8.082023423 21.771044
## 21   1.686299462  5.230567
## 22   1.931742422  7.918463
## 23   2.970551588  8.930210
## 24   2.696088971  9.913929
## 25   2.337057969  8.490298
## 26   3.211704686 10.246198
## 27   2.993485737  8.407554
## 28   3.141305094 12.239423
## 29   0.063535488  4.231198
## 30   1.195699383  5.914537
## 31   4.335679491 13.576897
## 32   1.013231457  5.903790
## 33   0.519115635  8.791728
## 34   5.193784233 16.371928
## 35   2.458051464  6.906298
## 36   1.581363518  6.112369
## 37   1.258560156  3.685271
## 38   2.509282006 11.132566
## 39   1.177369284  8.015315
## 40   2.258580068 12.294840
## 41   0.840729605  4.506440
## 42  14.422015152 37.346980
## 43   1.691443930  9.131246
## 44   0.451084013  1.273165
## 45   2.200677635  8.956381
## 46   4.496611385 13.076787
## 47   2.727468599 10.468491
## 48   1.152783336  4.391937
## 49   5.450551700 15.133240
## 50   2.624326085 15.412743
## 51   0.181182700  2.358367
## 52   0.612407700  3.299248
## 53   2.134426139  7.162627
## 54   0.627032512  2.892059
## 55   1.949280301  5.700664
## 56   3.775646631 14.098640
## 57   1.129177205  6.787550
## 58   5.153922658 14.847712
## 59   2.095391496  8.044475
## 60   2.048882684  6.714640
## 61   2.055739231  5.914148
## 62   0.569336921  9.945400
## 63   3.126103777  8.607377
## 64   0.084176586  7.601834
## 65   0.197261780  1.011401
## 66   0.197138624  3.461170
## 67   0.560770323  2.590245
## 68   0.591573917  2.120138
## 69   1.944859298  6.223975
## 70   1.848054448  6.494279
## 71   3.284809053 13.047825
## 72   3.239825535 10.272629
## 73   5.072291080 16.473134
## 74   3.043099247  8.733810
## 75   0.760028406  5.569288
## 76   0.477025934  4.935946
## 77   0.932974847  7.996663
## 78   0.084542903  1.224870
## 79   0.639353799  9.587242
## 80   1.287221841  6.599898
## 81   1.145126207  3.143113
## 82   0.432111005  1.342696
## 83   8.997346796 31.308371
## 84   3.714181106 14.771493
## 85   1.370917276  4.932154
## 86   2.878905252  7.384271
## 87   3.462307967  9.816682
## 88   2.489566600  8.213647
## 89   2.926601120 13.981065
## 90   3.074787957  9.580373
## 91   0.009198253  2.831010
## 92   2.217530936  6.726106
## 93   0.599940635  2.600885
## 94   2.384006014 10.739803
## 95   2.229857407  7.907911
## 96   0.134751782 11.420304
## 97   0.961337442  5.016904
## 98   3.140908678  9.787780
## 99   0.519892214  2.316506
## 100  3.713844575  9.812985

Ya teniendo las variables a estudiar, vamos a ver un gráfico de estas para ver si tienen o no distribución normal

plot(vector1, vector2)

Podemos ver que es poco probable que los datos tengan una distribución normal, pero vamos a probarlo y de paso verémos si tienen algún tipo de correlación.

H_0 = No hay correlación entre Vector 1 y Vector 2

H_A = Hay correlación entre Vector 1 y Vector 2

ks.test(vector1, "pnorm")
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  vector1
## D = 0.53816, p-value < 2.2e-16
## alternative hypothesis: two-sided
ks.test(vector2, "pnorm")
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  vector2
## D = 0.933, p-value < 2.2e-16
## alternative hypothesis: two-sided
cor.test(vector1, vector2, method = "spearman") #ya que ambas variables no presentan dN
## 
##  Spearman's rank correlation rho
## 
## data:  vector1 and vector2
## S = 33340, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##     rho 
## 0.79994

Como podemos observar, los datos efectivamente no presentan una distribución normal ya que el p-valor de ambas variables es menor a 0.05. Por otro lado, el p valor de la correlación es < 0.05, por lo tanto rechazamos H_0 y tomamos como cierta a H_A, en otras palabras, hay correlación entre Vector 1 y Vector 2. Además, tiene una correlación positiva.

Ahora, harémos un modelo de regresión de Theil-sen para intentar predecir el comportamiento de Vector 2 respecto a Vector 1. (install.packages(“mblm”))

library(mblm)

modelo_Ts <- mblm(vector2~vector1)

summary(modelo_Ts)
## 
## Call:
## mblm(formula = vector2 ~ vector1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2904 -0.9471  0.0033  1.5946  8.6177 
## 
## Coefficients:
##             Estimate    MAD V value Pr(>|V|)    
## (Intercept)   2.4785 1.2667    5041   <2e-16 ***
## vector1       2.4054 0.4955    5025   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.356 on 98 degrees of freedom

Donde Estimate esel valor esperado cuando vector 1 = 0 y, por cada aumento de una unidad de vector 1, vector 2 aumenta en 2.4054. En otras palabras, nuestro modelo sería: 2.4054x + 2.4785

Para mayor visualización, veamoslo en un gráfico:

plot(vector1, vector2, main = "Regresión de Theil-Sen")
abline(modelo_Ts, col = "red", lwd = 2)