Datos bivariados
Felipe González
Transformación para enderezar relaciones
Cuando consideramos la relación entre dos variables cuantitativas, muchas veces es posible transformar las variables para obtener una relación más simple y fácil de describir.
Ejemplo: desemepeño de CPUs
Utilizamos una base de datos de desemepeño de CPUs incluída en el paquete MASS, nos concentramos en dos variables: cuanto tarda un ciclo del procesador en nanosegundos (\(syct\)) y un índice de desempeño (\(estperf\)).
library(MASS)
ggplot(cpus, aes(x= syct, y = estperf)) +
geom_point(size = 1.2) +
geom_smooth(method = "loess", family = "symmetric", span = 0.4)
La relación entre estas dos variables es claramente inversa. Veamos los residuales para explorar si persiste una relación entre los residuales y la variable \(syct\) y para ver si el supuesto de homogeneidad es razonable en estos datos.
ajuste_loess_cpu <- loess(perf ~ syct, cpus, family = "symmetric", span = 0.2)
cpus$residual <- ajuste_loess_cpu$residual
cpus$ajustado <- ajuste_loess_cpu$fitted
residual_syct <- ggplot(cpus, aes(x = syct, y = residual)) +
geom_point(size = 1.2) +
geom_smooth(method = "loess", family = "symmetric", span = 1)
residual_ajuste <- ggplot(cpus, aes(x = ajustado, y = sqrt(abs(residual)))) +
geom_point(size = 1.2) +
geom_smooth(method = "loess", family = "symmetric", span = 1)
grid.arrange(residual_syct, residual_ajuste, ncol = 2)
Notamos que los residuales son considerablemente heterogéneos. Intentamos con una transformación: aplicamos logaritmo a la respuesta (perf).
ggplot(cpus, aes(x= syct, y = log(estperf))) +
geom_point(size = 1.2) +
geom_smooth(method = "loess", family = "symmetric", span = 0.4)
Y repetimos el análisis de residuales.
ajuste_loess_log_cpu <- loess(log(perf) ~ syct, cpus, family = "symmetric",
span = 0.2)
cpus$residual_log <- ajuste_loess_log_cpu$residual
cpus$ajustado_log <- ajuste_loess_log_cpu$fitted
residual_syct_log <- ggplot(cpus, aes(x = syct, y = residual_log)) +
geom_point(size = 1.2)+
geom_smooth(method = "loess", family = "symmetric", span = 1)
residual_ajuste_log <- ggplot(cpus, aes(x = ajustado_log,
y = sqrt(abs(residual_log)))) +
geom_point(size = 1.2) +
geom_smooth(method = "loess", family = "symmetric", span = 1)
grid.arrange(residual_syct_log, residual_ajuste_log, ncol = 2)
Notamos nuevamente que la relación entre syct y perf es inversa, quizá podemos transformar el factor syct en un intento de enderezar y simplificar esta relación, con el fin de obtener un modelo más simple. Usamos el inverso de cyst. Observemos que ya corregimos en buena parte la heterogeneidad de los residuales. Esta transformación, en cambio, tiene el propósito de enderezar la relación entre \(y\) y \(x\).
cpus$m_cycles_sec <- 1000 * (1 / cpus$syct) # millones de ciclos por segundo
ggplot(cpus, aes(x = m_cycles_sec, y = log(perf))) +
geom_point(size = 1.2) +
geom_smooth(method = "loess", family = "symmetric", degree = 1, span = 0.7)
Nótese que el inverso de syct es simplemente el número de ciclos por nanosegundo. Vemos también que la nueva relación es practicamente lineal.
En nuestro caso, tomando el inverso de nanosegundos por ciclo simplificamos la dependencia desempeño con este factor. Podemos revisar homogeneidad de residuales (obtenemos una gráfica muy similar a la que obtuvimos antes de transformar el factor).
ajuste_loess_trans <- loess(log(perf) ~ m_cycles_sec, cpus, family = "symmetric",
span = 0.2)
cpus$residual_trans <- ajuste_loess_trans$residual
cpus$ajustado_trans <- ajuste_loess_trans$fitted
residual_syct_trans <- ggplot(cpus, aes(x = m_cycles_sec, y = residual_trans)) +
geom_point(size = 1.2) +
geom_smooth(method = "loess", family = "symmetric", span = 1)
residual_ajuste_trans <- ggplot(cpus, aes(x = ajustado_trans,
y = sqrt(abs(residual_trans)))) +
geom_point(size = 1.2) +
geom_smooth(method = "loess", family = "symmetric", span = 1, degree = 1,
se = FALSE)
grid.arrange(residual_syct_trans, residual_ajuste_trans, ncol = 2)
Observamos que la dispersión es más homogénea (aunque para los cpus de peor desempeño la dispersión es un poco más baja), vemos también una consecuencia deseable: la relación entre respuesta y factor puede ahora describirse razonablemente mediante una línea recta.
Series de tiempo
Podemos usar el suavizamiento loess para entender y describir el comportamiento de series de tiempo, en las cuales intentamos entender la dependencia de una serie de mediciones indexadas por el tiempo. Típicamente es necesario utilizar distintas componentes para describir exitosamente una serie de tiempo, y para esto usamos distintos tipos de suavizamientos. Veremos que distintas componentes varían en distintas escalas de tiempo (unas muy lentas, cono la tendencia, otras más rapidamente, como variación quincenal, etc.).
En el siguiente ejemplo consideramos la ventas semanales de un producto a lo largo de 5 años. Veamos que existe una tendencia a largo plazo (crecimientos anuales) y también que existen patrones de variación estacionales.
ventas <- read.csv("./datos/ventas_semanal.csv")
head(ventas) period sales.kg
1 1 685.537
2 2 768.234
3 3 894.643
4 4 773.501
5 5 954.600
6 6 852.853ggplot(ventas, aes(x = period, y = sales.kg)) + geom_line(size = 0.3)
Intentaremos usar suavizamiento para capturar los distintos tipos de variación que observamos en la serie. En primer lugar, si suavizamos poco (por ejemplo \(\aplha = 0.1\)), vemos que capturamos en parte la tendencia y en parte la variación estacional.
ggplot(ventas, aes(x = period, y = log(sales.kg))) +
geom_line(size = 0.3) +
geom_smooth(method = "loess", span = 0.1, degree = 1, se = FALSE, size = 1,
color = "red")
Es mejor comenzar capturando la tendencia, y poco de la componente estacional:
ggplot(ventas, aes(x = period, y = log(sales.kg))) +
geom_line(size = 0.3) +
geom_smooth(method = "loess", span = 0.3, degree = 1, se = FALSE, size = 1,
color = "red")
ajuste.trend.1 <- loess(log(sales.kg) ~ period, ventas, span = 0.3, degree = 1)
ventas$trend.1 <- ajuste.trend.1$fitted
ventas$res.trend.1 <- ajuste.trend.1$residualsAhora calculamos los residuales de este ajuste e intentamos describirlos mediante un suavizamiento más fino. Verificamos que hemos estimado la mayor parte de la tendencia, e intentamos capturar la variación estacional de los residuales.
ggplot(ventas, aes(x = period, y = res.trend.1)) +
geom_line(size = 0.3) +
geom_smooth(method = "loess", span = 0.15, degree = 1, se = FALSE, size = 1,
color = "red")
ajuste.est1.1 <- loess(res.trend.1 ~ period, ventas, span = 0.15, degree = 1)
ventas$est1.1 <- ajuste.est1.1$fitted
ventas$res.est1.1 <- ajuste.est1.1$residualsY graficamos los residuales obtenidos después de ajustar el componente estacional para estudiar la componente de mayor frecuencia.
ggplot(ventas, aes(x = period, y = res.est1.1)) +
geom_line(size = 0.3) +
geom_smooth(method = "loess", span = 0.06, degree = 1, se = FALSE, size = 1,
color = "red")
ajuste.est2.1 <- loess(res.est1.1 ~ period, ventas, span = 0.06, degree = 1)
ventas$est2.1 <- ajuste.est2.1$fitted
ventas$res.est2.1 <- ajuste.est2.1$residualsAhora que tenemos nuestra primera estimación de cada una de las componentes, podemos regresar a hacer una mejor estimación de la tendencia. La ventaja de volver es que ahora podemos suavizar más sin que en nuestra muestra compita tanto la variación estacional. Por tanto podemos suavizar menos:
ventas$sales.sin.est.1 <- log(ventas$sales.kg) - ajuste.est1.1$fitted -
ajuste.est2.1$fitted
ggplot(ventas, aes(x = period, y = sales.sin.est.1)) +
geom_line(size = 0.3) +
geom_smooth(method = "loess", span = 0.08, degree = 1, se = FALSE, size = 1,
color = "red")
ajuste.trend.2 <- loess(sales.sin.est.1 ~ period, ventas, span = 0.08, degree = 1)
ventas$trend.2 <- ajuste.trend.2$fitted
ventas$res.trend.2 <- log(ventas$sales.kg) - ventas$trend.2Y ahora nos concentramos en la componente anual.
ventas$sales.sin.est.2 <- log(ventas$sales.kg) - ajuste.trend.2$fitted -
ajuste.est2.1$fitted
ggplot(ventas, aes(x = period, y = sales.sin.est.2)) +
geom_line(size = 0.3) +
geom_smooth(method = "loess", span = 0.2, degree = 2, se = FALSE, size = 1,
color = "red")
ajuste.est1.2 <- loess(sales.sin.est.2 ~ period, ventas, span = 0.15, degree = 1)
ventas$est1.2 <- ajuste.est1.2$fitted
ventas$res.est1.2 <- ajuste.est1.2$residualsFinalmente volvemos a ajustar la componente de frecuencia más alta:
ventas$sales.sin.est.3 <- log(ventas$sales.kg) - ajuste.trend.2$fitted -
ajuste.est1.2$fitted
ggplot(ventas, aes(x = period, y = sales.sin.est.3)) +
geom_line(size = 0.3) +
geom_smooth(method = "loess", span = 0.06, degree = 1, se = FALSE, size = 1,
color = "red")
ajuste.est2.2 <- loess(sales.sin.est.3 ~ period, ventas, span = 0.06, degree = 1)
ventas$est2.2 <- ajuste.est2.2$fitted
ventas$res.est2.2 <- ajuste.est2.2$residualsVerificamos nuestra descomposición y visualizamos el ajuste:
ventas$log.sales <- log(ventas$sales.kg)
ventas.2 <- dplyr::select(ventas, period, trend.2, est1.2, est2.2, res.est2.2,
log.sales)
max(abs(apply(ventas.2[, 2:4], 1, sum) - ventas.2$log.sales))[1] 0.240316ventas.2.m <- gather(ventas.2, componente, valor, -period)
ventas.2.m.c <- ventas.2.m %>%
group_by(componente) %>%
mutate(
valor.c = valor - mean(valor)
)
ggplot(ventas.2.m.c, aes(x = period, y = valor.c)) +
geom_vline(x = c(0, 52 - 1, 52 * 2 - 1, 52 * 3 - 1, 52 * 4 - 1), color = "gray") +
geom_line(size = 0.3) +
facet_wrap(~ componente, ncol = 1)
Y vemos que es razonable describir los residuales con una distribución normal (con desviación estándar alrededor de 8% sobre el valor ajustado):
sd(ventas$res.est2.2)[1] 0.08093687ventas.ord <- arrange(ventas, res.est2.2)
ventas.ord$q.normal <- qnorm((1:nrow(ventas) - 0.5) / nrow(ventas))
ggplot(ventas.ord, aes(x = q.normal, y = res.est2.2)) +
geom_point(size = 1.2) +
geom_smooth(method = "lm")
Aún no estamos convencidos de que podemos hacer pooling de los residuales. Para checar esto, vemos la gráfica de dependencia de residuales.
ggplot(ventas, aes(x = period, y = res.est2.2)) + geom_line(size = 0.3) +
geom_point(size = 1.2)
¿Puedes notar alguna dependencia en los residuales?
library(nullabor)
ventas.res <- dplyr::select(ventas, period, res.est2.2)
ventas.null <- lineup(null_dist(var = 'res.est2.2', dist = 'normal',
params = list(mean = 0, sd = 0.08)), n = 10, ventas.res)decrypt("OlCE bQTQ Aw GWPATAWw 6")ggplot(ventas.null, aes(x = period, y = res.est2.2)) +
facet_wrap(~ .sample, ncol = 2) +
geom_line(size = 0.3) +
geom_vline(x = c(0, 52 - 1, 52 * 2 - 1, 52 * 3 - 1, 52 * 4 - 1), color = "gray") +
geom_point(size = 1.2) 
Hay dos cosas que nos falta explicar, en primer lugar, las caídas alrededor de principios/finales de cada año (que son de hasta -0.2), y segundo que esta gráfica parece oscilar demasiado. La estructura que aún no hemos explicado se debe a que las semanas que caen en quincena tienden a tener compras más grandes que las que están justo antes de quincena o fin de mes.
Por el momento detendremos el análisis aquí y explicamos un proceso iterativo para proceder en nuestro análisis exploratorio:
Iterando ajuste de loess. Cuando queremos ajustar con tres componentes: tendencia, estacionalidad y residuales, podemos seguir el siguiente proceso,
Ajustar la primera componente a los datos (tendencia).
Ajustar la segunda componente a los residuales del paso anterior (estacionalidad).
Restar de los datos originales la segunda componente ajustada (estacionalidad).
Ajustar a los residuales del paso anterior una nueva componente (tendencia).
Restar a los datos originales la componente ajustada en el paso anterior.
Ajustar a los residuales del paso anterior una nueva componente (estacionalidad).
- Checar ajuste y si es necesario iterar de 3 a 6 con las nuevas componentes.
La idea es que cada componente compite para explicar los datos (cada una gana más al bajar el parámetro \(\alpha\)). El conflicto es que si suavizamos mucho cada componente (por ejemplo la tendencia), entonces parte de la variación que debería ir en ella queda en los residuales, y se intenta ajustar posteriormente por una componente distinta (estacionalidad). Sin embargo, si suavizamos poco, entonces parte de la variación de la segunda componente es explicada por el ajuste de la primera. Entonces, la solución es ir poco a poco adjudicando variación a cada componente. En nuestro ejemplo de arriba, podemos comenzar suavizando de menos el primer ajsute de la tendencia, luego ajustar estacionalidad, restar a los datos originales esta estacionalidad, y ajustar a estos datos una componente más suave de tandencia. Es posible suavizar más la tendencia justamente porque ya hemos eliminado una buena parte de la estacionalidad.
Ahora, si vemos como se comportan los residuales según el día donde comienza la semana, vemos el patrón que explicamos antes:
dat.tot <- read.csv(file = "./datos/cereal_tot.csv")
head(dat.tot) period sales.kg year month day week quincena price wd wpd
1 1 685.537 2002 1 6 1 0 25.83813 5529.553 959.7047
2 2 768.234 2002 1 13 2 0 26.03183 5558.863 932.0702
3 3 894.643 2002 1 20 3 1 25.64726 5640.989 1629.9504
4 4 773.501 2002 1 27 4 0 25.58614 5665.832 1580.0875
5 5 954.600 2002 2 3 5 1 25.49396 5701.529 1642.4942
6 6 852.853 2002 2 10 6 0 25.66043 5658.477 1243.6257
Variants TRPs adstocks trend1 trend2 y2002
1 85 0.00 6000.000 1 0.6931472 85
2 87 472.11 5872.110 2 1.0986123 87
3 87 505.64 5790.539 3 1.3862944 87
4 87 1054.22 6265.705 4 1.6094379 87
5 87 548.09 6187.225 5 1.7917595 87
6 86 745.85 6314.352 6 1.9459101 86ventas.day <- inner_join(ventas, dplyr::select(dat.tot, period, day))Joining by: "period"ggplot(ventas.day, aes(x = day, y = res.est2.2)) +
geom_point() +
ylab("residual") +
geom_smooth(method = "loess", span = 0.06, degree = 1, se = FALSE, size = 1,
color = "red")Warning in loop_apply(n, do.ply): pseudoinverse used at 0.85Warning in loop_apply(n, do.ply): neighborhood radius 1.15Warning in loop_apply(n, do.ply): reciprocal condition number 0Warning in loop_apply(n, do.ply): There are other near singularities as
well. 1
Podemos hacer un ajuste loess con estos residuales:
ajuste.quincenas <- loess(res.est2.2 ~ day, data = ventas.day)
ventas$quincena <- ajuste.quincenas$fitted
ventas$res.final <- ajuste.quincenas$residuals
sd(ventas$quincena)[1] 0.04302381sd(ventas$res.final)[1] 0.06359823ggplot(ventas, aes(x = period, y = res.final)) +
geom_line(size = 0.3) +
geom_point(size = 1.2)
ventas.2 <- dplyr::select(ventas, period, trend.2, est1.2, est2.2, quincena,
res.final, log.sales)
max(abs(apply(ventas.2[, 2:4], 1, sum) - ventas.2$log.sales))[1] 0.240316ventas.2.m <- gather(ventas.2, componente, valor, -period)
ventas.2.m.c <- ventas.2.m %>%
group_by(componente) %>%
mutate(
valor.c = valor - mean(valor)
)
ggplot(ventas.2.m.c, aes(x = period, y = valor.c)) +
geom_vline(x = c(0, 52 - 1, 52 * 2 - 1, 52 * 3 - 1, 52 * 4 - 1), color = "gray") +
geom_line(size = 0.3) +
facet_wrap(~ componente, ncol = 1)