Tema: Trabajo corto 6

Curso: Análisis de regresión

Profesor: Jesús Eduardo Gamboa

UNIVERSIDAD NACIONAL AGRARIA LA MOLINA

1. INTRODUCCIÓN

CONTEXTO DE LOS DATOS

¿Qué son los animes?

Es un estilo de animación tradicional o por computadora de origen japones el cual cuenta con una gran variedad de generos y tramas. Este producto es consumido tanto en Japon como alrededor del mundo y ha venido teniendo un gran impacto en los hispanohablantes.

Segun la página de noticias ReporteASIA, comenta que En 2020 el mercado del anime (dibujo animado japonés) a nivel global fue valuado en US 24.23 mil millones de dolares y se prevé que alcance los US 43.73 mil millones de dolares en 2027 con una tasa compuesta anual del 8,8%. Dicha información fue extraida del AJA (Association of Japanese Animation).

Se trabajo con la base de datos de una página bastante reconocida de animes y mangas, llamada “My Anime List” la cual fue publicado en keaggle el 3 de julio del 2022. Para su recolección de los datos tan solo fue necesario hacer uso del API de la página de animes.

La extracción de datos fue realizada por Andreu Vall Hernàndez, estudiante de matemática en la Universidad de Barcelona.

FUENTES

Para mas informacion, revisar los siguientes hipervinculos.

ReporteASIA

Keaggle

My Anime List

PREPARACIÓN DE DATOS

Lectura y limpieza de datos

Empezamos leyendo los datos a trabajar y eliminando todos los valores nulos (Lo ultimo solo por cuestión de practidad)

Tambien nos quedamos tan solo con 4 variables de las 31.

Las cuales son:

Nombre Tipo de dato Unidades Rango (MXN) Descripción
titulo Nintendo cualitativo nominal - Es el nombre del anime.
puntaje cuantitativo continuo - 0 a 10 Es el puntaje dado por los usuarios de MAL.
cantidad_episodios Cuantitativo discreto episodio 0 a infinito Es la cantidad de episodios que tiene un anime.
popularidad Cuantitativo discreto vistas 0 a infinito Es la cantidad de personas que han visto el anime. 
anime = read.csv("anime.csv") |> 
  drop_na() |> 
  select("title","mean","num_episodes","popularity") |> 
  rename(titulo = title, puntaje = mean, cantidad_episodios = num_episodes, popularidad = popularity)

Muestreo

Dado que originalmente hay 20602 observaciones, optamos por realizar un muestreo simple para extraer tan solo 30 datos.

N<-nrow(anime) 
n<-30
RNGkind(sample.kind="Rounding")
set.seed(20)
m1<-S.SI(N,n)
mues<-anime[m1,]

Visualización de la tabla de datos

En la siguiente tabla se puede visualizar los valores obtenidos en la muestra de tamaño 30.

mues |> datatable(option=list(scrollX=T))

2. DESARROLLO DE PREGUNTAS

1. Utilizar un conjunto de datos con 30 observaciones y 2 variables explicativas (y, x1, x2)

Selección de variables

Su seleccion fue solo por cuestion de practicida, no hubo criterios previos para escoger alguno u orden de variables.

x1 = mues$puntaje #explicativa
x2 = mues$popularidad #explicativa
y = mues$cantidad_episodios #predictiva

Formación modelo

\[\widehat{Y} = \beta_0+X\beta_1+X\beta2\]

lm(y ~ x1+x2, data = mues) -> modelo

2. Utilizar los residuales estandarizados o estudentizados, según sea lo adecuado, para identificar posibles outliers.

RESIDUALES ESTANDARIZADOS

Como no sabemos si el residual es grande o pequeño porque sus unidades son relativas, entonces aplicamos el residual estadarizados, el cual no tiene unidades. Lo cual nos permite tener una escala.

\[d_i=\frac{e_i}{\widehat{\sigma}\sqrt{1-h_i}}\] donde: - \(e_i\) es el i-ésimo residual - \(h_ii\) es el i-ésimo leverage (i-ésimo termino de la matriz hat) - \(\widehat{\sigma}^2\) es la estimación de la desviación estándar, recordando que para un modelo con n observaciones y k coeficientes de regresión:

\[\widehat{\sigma}^2=\frac{1}{n-k}\sum^{n}_{j=1}e^2_j=\widehat{CME}\]

Tener en cuenta que: - La media de un residual estandarizado es 0 y su varianza es 1. Con lo que se puede asumir normalidad. - Valores por encima de 2, o por debajo de -2 son considerados outliers.

Obteniendo el resultado de forma operativa

modelo |> resid() -> r
summary(modelo)$sigma -> s
modelo |> model.matrix() -> X
X %*% solve(t(X) %*% X) %*% t(X) -> H
res_stand <- (r/(s*sqrt(1-diag(H)))) 
res_stand |> round(3)
##     13    208    764   1032   1138   1162   1451   2692   3783   4661   5038 
##  0.142  0.083 -0.634 -0.766 -0.738  0.326  0.095 -0.370  0.252 -0.097 -0.220 
##   5039   5054   5146   6520   6551   6599   7213   7424   8533   8721   8948 
##  0.179 -0.363  0.313 -0.664 -0.375 -0.058  2.983  3.581 -0.410 -1.027 -0.778 
##   9359   9675  10347  10607  10621  10666  11402  11926 
##  0.362 -0.525  0.124 -0.803  0.819 -0.383 -0.485 -0.701

Comprobando con la formula directa

modelo |> rstandard() |> round(3) |> abs()
##    13   208   764  1032  1138  1162  1451  2692  3783  4661  5038  5039  5054 
## 0.142 0.083 0.634 0.766 0.738 0.326 0.095 0.370 0.252 0.097 0.220 0.179 0.363 
##  5146  6520  6551  6599  7213  7424  8533  8721  8948  9359  9675 10347 10607 
## 0.313 0.664 0.375 0.058 2.983 3.581 0.410 1.027 0.778 0.362 0.525 0.124 0.803 
## 10621 10666 11402 11926 
## 0.819 0.383 0.485 0.701
(modelo |> rstandard() |> round(3) |> abs()) > 2
##    13   208   764  1032  1138  1162  1451  2692  3783  4661  5038  5039  5054 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##  5146  6520  6551  6599  7213  7424  8533  8721  8948  9359  9675 10347 10607 
## FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
## 10621 10666 11402 11926 
## FALSE FALSE FALSE FALSE

Gráficando

Podemos observar que tanto la observacion 18 y 19 son outliers.

modelo |> ols_plot_resid_stand()

Residuales estudentizados

Dado que hemos encontrado dos outiliers, es mejor trabajar con los residuales estudentizados (externamente).

Es bastante similar al estandarizado o tambien llamado estudentizado internamente, solo que se quita e i-ésimo termino en los residuales \[t_i=\frac{e_i}{\widehat{\sigma}_{(i)}\sqrt{1-h_{ii}}}=d_i(\frac{n-k-1}{n-k-d^2_i})^{0.5}\] \[\widehat{\sigma}^2_{(i)}=\frac{1}{n-k-1}\sum^{n}_{j=1,j\neq 1}e^2_i=\widehat{CME_{(i)}}\]

Obteniendo el resultado de forma operativa

lm(y[-18]~x1[-18]+x2[-18]) -> modelo_1
summary(modelo_1)$sigma -> s_1
r[18]/(s_1*sqrt(1-H[18,18]))
##     7213 
## 3.574184
lm(y[-19]~x1[-19]+x2[-19]) -> modelo_2
summary(modelo_2)$sigma -> s_2
r[19]/(s_2*sqrt(1-H[19,19]))
##     7424 
## 4.850034
length(x1) -> n; 3 -> k
(res_stand*sqrt((n-k-1)/(n-k-res_stand**2))) |> round(3)
##     13    208    764   1032   1138   1162   1451   2692   3783   4661   5038 
##  0.140  0.081 -0.627 -0.760 -0.732  0.321  0.094 -0.364  0.248 -0.095 -0.216 
##   5039   5054   5146   6520   6551   6599   7213   7424   8533   8721   8948 
##  0.175 -0.357  0.307 -0.657 -0.368 -0.057  3.574  4.850 -0.403 -1.028 -0.772 
##   9359   9675  10347  10607  10621  10666  11402  11926 
##  0.356 -0.518  0.122 -0.797  0.814 -0.377 -0.478 -0.694

Comprobando con la formula directa

library(MASS)

modelo |> studres() |> round(3) |> abs()
##    13   208   764  1032  1138  1162  1451  2692  3783  4661  5038  5039  5054 
## 0.140 0.081 0.627 0.760 0.732 0.321 0.094 0.364 0.248 0.095 0.216 0.175 0.357 
##  5146  6520  6551  6599  7213  7424  8533  8721  8948  9359  9675 10347 10607 
## 0.307 0.657 0.368 0.057 3.574 4.850 0.403 1.028 0.772 0.356 0.518 0.122 0.797 
## 10621 10666 11402 11926 
## 0.814 0.377 0.478 0.694
(modelo |> studres() |> round(3) |> abs()) > 2
##    13   208   764  1032  1138  1162  1451  2692  3783  4661  5038  5039  5054 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##  5146  6520  6551  6599  7213  7424  8533  8721  8948  9359  9675 10347 10607 
## FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
## 10621 10666 11402 11926 
## FALSE FALSE FALSE FALSE

Gráficando

modelo |> ols_plot_resid_stud()

3. Identificar las posibles observaciones leverage

modelo %>% model.matrix -> X
X %*% solve(t(X) %*% X) %*% t(X) -> H
#H

n <- length(x1)
k <- 3
diag(H)
##         13        208        764       1032       1138       1162       1451 
## 0.31948632 0.10992568 0.09434551 0.10474281 0.09801994 0.15628162 0.08260830 
##       2692       3783       4661       5038       5039       5054       5146 
## 0.10291716 0.07466751 0.15291144 0.04467978 0.04743484 0.04776001 0.08390564 
##       6520       6551       6599       7213       7424       8533       8721 
## 0.03881203 0.03473863 0.03775838 0.08036914 0.03668708 0.06741522 0.16100909 
##       8948       9359       9675      10347      10607      10621      10666 
## 0.07849395 0.06514705 0.06440088 0.23523912 0.11400120 0.08569817 0.08441107 
##      11402      11926 
## 0.12301523 0.17311721
2*k/n
## [1] 0.2
diag(H) > 2*k/n
##    13   208   764  1032  1138  1162  1451  2692  3783  4661  5038  5039  5054 
##  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
##  5146  6520  6551  6599  7213  7424  8533  8721  8948  9359  9675 10347 10607 
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE 
## 10621 10666 11402 11926 
## FALSE FALSE FALSE FALSE
data.frame(leverage = diag(H)>2*k/n,
           outlier = abs(modelo |> rstudent())>2)
##       leverage outlier
## 13        TRUE   FALSE
## 208      FALSE   FALSE
## 764      FALSE   FALSE
## 1032     FALSE   FALSE
## 1138     FALSE   FALSE
## 1162     FALSE   FALSE
## 1451     FALSE   FALSE
## 2692     FALSE   FALSE
## 3783     FALSE   FALSE
## 4661     FALSE   FALSE
## 5038     FALSE   FALSE
## 5039     FALSE   FALSE
## 5054     FALSE   FALSE
## 5146     FALSE   FALSE
## 6520     FALSE   FALSE
## 6551     FALSE   FALSE
## 6599     FALSE   FALSE
## 7213     FALSE    TRUE
## 7424     FALSE    TRUE
## 8533     FALSE   FALSE
## 8721     FALSE   FALSE
## 8948     FALSE   FALSE
## 9359     FALSE   FALSE
## 9675     FALSE   FALSE
## 10347     TRUE   FALSE
## 10607    FALSE   FALSE
## 10621    FALSE   FALSE
## 10666    FALSE   FALSE
## 11402    FALSE   FALSE
## 11926    FALSE   FALSE
modelo |> plot(which=5)

modelo |> ols_plot_resid_lev()

4. Identificar los posibles valores influenciales según distancia de Cook (interpretar para al menos 2 alfa)

aplicacndo la distancia de cook

modelo |> cooks.distance()
##           13          208          764         1032         1138         1162 
## 3.170415e-03 2.818784e-04 1.396655e-02 2.289653e-02 1.974457e-02 6.572109e-03 
##         1451         2692         3783         4661         5038         5039 
## 2.731141e-04 5.249088e-03 1.711508e-03 5.636401e-04 7.545250e-04 5.289251e-04 
##         5054         5146         6520         6551         6599         7213 
## 2.204134e-03 2.982389e-03 5.929098e-03 1.682676e-03 4.370543e-05 2.591329e-01 
##         7424         8533         8721         8948         9359         9675 
## 1.628067e-01 4.043441e-03 6.744094e-02 1.718341e-02 3.042199e-03 6.330607e-03 
##        10347        10607        10621        10666        11402        11926 
## 1.585371e-03 2.763656e-02 2.094643e-02 4.503230e-03 1.100374e-02 3.431471e-02

gráficando

modelo |> ols_plot_cooksd_bar()

interpretación de los alfas

alfas = c(0.01,0.05,0.1,0.25,0.5,0.75,0.90,0.95,0.99)
data.frame(alfas,valorF=qf(alfas,k,n-k))
##   alfas     valorF
## 1  0.01 0.03767053
## 2  0.05 0.11592389
## 3  0.10 0.19335641
## 4  0.25 0.40584945
## 5  0.50 0.80894686
## 6  0.75 1.45096086
## 7  0.90 2.29871191
## 8  0.95 2.96035132
## 9  0.99 4.60090689

  • Si \(D_i \geqslant 0.03767053 \rightarrow \widehat{\beta}\) se mueve al menos 1% hacia el borde de la región de confianza

  • Si \(D_i \geqslant 0.80894686 \rightarrow \widehat{\beta}\) se mueve al menos 50% hacia el borde de la región de confianza

  • Si \(D_i \geqslant 2.96035132 \rightarrow \widehat{\beta}\) se mueve al menos 95% hacia el borde de la región de confianza

5. Identificar los posibles valores influenciales según DFBETAS

Aplicando dfbetas

modelo |> dfbetas()
##        (Intercept)           x1           x2
## 13    -0.069073882  0.077685289  0.018076410
## 208    0.009869707 -0.006018465 -0.021063653
## 764   -0.016958063 -0.011441717  0.111998923
## 1032   0.066217256 -0.100728365  0.071803328
## 1138   0.047987254 -0.080558782  0.077369132
## 1162   0.105606821 -0.092156029 -0.121710308
## 1451   0.006988634 -0.003128648 -0.017947283
## 2692  -0.089977767  0.077666419  0.100660950
## 3783   0.048344960 -0.041355921 -0.051636055
## 4661  -0.036715175  0.034553937  0.029925007
## 5038  -0.022724393  0.018004552  0.023353292
## 5039   0.021202766 -0.017395233 -0.020846180
## 5054  -0.044000963  0.036286337  0.042812622
## 5146  -0.061673780  0.069270315  0.031787510
## 6520   0.039961519 -0.049152386 -0.037975484
## 6551  -0.018029701  0.013371051  0.005825313
## 6599   0.003001242 -0.003760780 -0.003136895
## 7213  -0.749298448  0.789454565  0.661220457
## 7424  -0.008123085  0.045354857  0.239008885
## 8533  -0.060999282  0.061185591  0.006742995
## 8721   0.373387336 -0.377004341 -0.355303808
## 8948   0.132829463 -0.132947924 -0.169012538
## 9359  -0.033926681  0.032801701  0.064110268
## 9675  -0.020091623  0.023446504 -0.051610210
## 10347  0.053820089 -0.055755563 -0.016479053
## 10607  0.159213609 -0.150752850 -0.240077512
## 10621  0.036470991 -0.047002821  0.107612178
## 10666 -0.005059619  0.009931256 -0.058493274
## 11402 -0.032955649  0.043311817 -0.078890734
## 11926  0.083541341 -0.062374758 -0.247613935
(abs(modelo |> dfbetas())) >= 2/sqrt(n)
##       (Intercept)    x1    x2
## 13          FALSE FALSE FALSE
## 208         FALSE FALSE FALSE
## 764         FALSE FALSE FALSE
## 1032        FALSE FALSE FALSE
## 1138        FALSE FALSE FALSE
## 1162        FALSE FALSE FALSE
## 1451        FALSE FALSE FALSE
## 2692        FALSE FALSE FALSE
## 3783        FALSE FALSE FALSE
## 4661        FALSE FALSE FALSE
## 5038        FALSE FALSE FALSE
## 5039        FALSE FALSE FALSE
## 5054        FALSE FALSE FALSE
## 5146        FALSE FALSE FALSE
## 6520        FALSE FALSE FALSE
## 6551        FALSE FALSE FALSE
## 6599        FALSE FALSE FALSE
## 7213         TRUE  TRUE  TRUE
## 7424        FALSE FALSE FALSE
## 8533        FALSE FALSE FALSE
## 8721         TRUE  TRUE FALSE
## 8948        FALSE FALSE FALSE
## 9359        FALSE FALSE FALSE
## 9675        FALSE FALSE FALSE
## 10347       FALSE FALSE FALSE
## 10607       FALSE FALSE FALSE
## 10621       FALSE FALSE FALSE
## 10666       FALSE FALSE FALSE
## 11402       FALSE FALSE FALSE
## 11926       FALSE FALSE FALSE

Gráficando

modelo |> ols_plot_dfbetas()