Evolución de la eficiencia aerodinámica de los automóviles argentinos 1960-2000. Descartando el resto de los factores lo único que limita la velocidad máxima dada cierta potencia es la resistencia del viento.

library(tidyverse)
library(ggplot2)
library(plotly)
library(akima)

Levanto la info del archivo. Verán que sigo con R 3.6.

df_autos_4 <- read.csv("df_autos_4.csv", stringsAsFactors = FALSE)
head(df_autos_4)
##   listavel listaanio listapotencia
## 1      131      1961            40
## 2       86      1963            34
## 3       85      1964            14
## 4       92      1966            14
## 5      114      1972            36
## 6      112      1971            35

Voy a armar seis categorías de potencia: 50, 70 hasta 150hp y dividir todo por décadas. Saco el promedio entre cada uno y después lo llevo (por regla de tres simple) hacia la definición de categoría. Ej: si para la categoría de 50hp de 1960 el promedio de los autos de 40 a 60hp es de 48hp y 123km/h se ajusta a 123*(50/48)=128,1 km/h. Con esto buscamos comparar entre un auto “ideal” de 1960 con 50hp exactos y eventualmente (por ejemplo) otro de 1980 con las mismas condiciones y características.

df_autos_clasifpot <- df_autos_4 %>% select(anio = listaanio, velmax = listavel, potencia = listapotencia) %>% drop_na() %>% filter(potencia >= 40 & potencia < 160) %>% mutate(
  clasifpotencia = case_when(
    potencia >= 40 & potencia < 60 ~ 50,
    potencia >= 60 & potencia < 80 ~ 70,
    potencia >= 80 & potencia < 100 ~ 90,
    potencia >= 100 & potencia < 120 ~ 110,
    potencia >= 120 & potencia < 140 ~ 130,
    potencia >= 140 & potencia < 160 ~ 150
   ),
  decada = (trunc(anio/10))*10
) %>% group_by(decada, clasifpotencia) %>% summarize(promvel = mean(velmax), prompot = mean(potencia)) %>%
  mutate(
    velajust = promvel*(clasifpotencia / prompot),
    relvelpot = promvel/prompot,
    potchar = as.character(sprintf("%03d",clasifpotencia)),
    decadachar = as.character(decada)
  )

Podemos ver que sale a simple vista.

ggplot() + 
  geom_line(data = df_autos_clasifpot, aes(x = potchar, y = velajust, group = decadachar, color = decadachar), size = 1.25 ) +
  labs(x = "Potencia (cv)",
       y = "Velocidad máxima (km/h)",
       colour = "Década"
       )

Es un gráfico con mucho ruido para una lectura cómoda. Veamos la relación entre velocidad y potencia:

ggplot() + 
  geom_line(data = df_autos_clasifpot, aes(x = potchar, y = relvelpot, group = decadachar, color = decadachar), size = 1.25 ) +
   labs(x = "Potencia (cv)",
       y = "Relación vel/pot kmh/cv",
       colour = "Década"
       )

Acá se ve mejor como se va dibujando la “curva” que todavía sólo existe como un puñado de líneas entre puntos y mucha imaginación. Con esta forma vemos que se parece mucho a un logaritmo, veamos si se lo puede ajustar. Volvemos a usar los datos originales sin la parte de corregir para aproximar a la definición de cada categoría. O sea por cada rango se toma el promedio y listo. Esta información se arma en pares de listas x/y (x1960-y1960…) y se utiliza para crear el modelo con lm.

Debe haber alguna forma fácil de hacer todo esto.

df_autos_1960 <- df_autos_clasifpot %>% filter(decada == 1960)
df_autos_1970 <- df_autos_clasifpot %>% filter(decada == 1970)
df_autos_1980 <- df_autos_clasifpot %>% filter(decada == 1980)
df_autos_1990 <- df_autos_clasifpot %>% filter(decada == 1990)
df_autos_2000 <- df_autos_clasifpot %>% filter(decada == 2000)

x1960 <- df_autos_1960$prompot
y1960 <- df_autos_1960$promvel
x1970 <- df_autos_1970$prompot
y1970 <- df_autos_1970$promvel
x1980 <- df_autos_1980$prompot
y1980 <- df_autos_1980$promvel
x1990 <- df_autos_1990$prompot
y1990 <- df_autos_1990$promvel
x2000 <- df_autos_2000$prompot
y2000 <- df_autos_2000$promvel

fit1960  <- lm(y1960~log(x1960))
fit1970  <- lm(y1970~log(x1970))
fit1980  <- lm(y1980~log(x1980))
fit1990  <- lm(y1990~log(x1990))
fit2000  <- lm(y2000~log(x2000))

xx <- seq(50,180, length=10)

Una vez armado el modelo lo dibujamos. En el primer plot hay que especificar que no le ponga etiquetas a los ejes por que si no después quedan superpuestos.

plot(x1960,y1960, col = "green", xlab = "", ylab = "", ylim = c(100,200))
points(x1970,y1970, col = "blue")
points(x1980,y1980, col = "grey")
points(x1990,y1990, col = "red")
points(x2000,y2000, col = "orange")
lines(xx, predict(fit1960, data.frame(x1960=xx)), col="green")
lines(xx, predict(fit1970, data.frame(x1970=xx)), col="blue")
lines(xx, predict(fit1980, data.frame(x1980=xx)), col="grey")
lines(xx, predict(fit1990, data.frame(x1990=xx)), col="red")
lines(xx, predict(fit2000, data.frame(x2000=xx)), col="orange")
title(xlab = "Potencia (hp)", ylab = "Velocidad máxima (km/h)")

Con estos datos se podría afirmar que hubo una evolución en términos de eficiencia aerodinámica de los automóviles a lo largo de las cinco décadas estudiadas. Un automóvil promedio con 80hp en las décadas del 60 y 70 estaba en aproximadamente 140km/h mientras que con la misma potencia uno de la década del 90 o 2000 podía hacer más de 160km/h. Sólo desentona un poco la década de 1980 pero es posible atribuir este fenómeno a que durante ese período se dio una transición en el parque automotor y los automóviles más grandes (ergo más potentes) como tienen ciclos de producto más largos hayan quedado rezagados en la tecnología mientras que comenzaron a producirse automóviles pequeños con tecnología más actualizada.

Hay un trade off entre la relación de los datos visualizados contra los datos reales (que son los del archivo) y un gráfico “cómodo” y simple de ver como este. Queda a criterio de cada analista elegir el mejor enfoque de acuerdo a la situación. En este caso estoy haciendo un modelo sobre conjuntos de promedios. Hay dos niveles de abstracción entre el dato y la presentación (primero el promedio y después el logaritmo). En el gráfico anterior había uno solo (el promedio).

¿Se pueden hacer modelos sobre todo el conjunto de datos? Por supuesto. Veamos por que no siempre es una buena idea. Notar que como acá no hay promedios el rango es directamente 50-150hp no 40-160 como el anterior.

df_autos_clasiftot <- df_autos_4 %>% select(anio = listaanio, velmax = listavel, potencia = listapotencia) %>% drop_na() %>% filter(potencia >= 50 & potencia < 150) %>% mutate(
  decada = (trunc(anio/10))*10
)

Insisto: tiene que haber alguna forma mejor de hacer esto.

df_autos_1960t <- df_autos_clasiftot %>% filter(decada == 1960)
df_autos_1970t <- df_autos_clasiftot %>% filter(decada == 1970)
df_autos_1980t <- df_autos_clasiftot %>% filter(decada == 1980)
df_autos_1990t <- df_autos_clasiftot %>% filter(decada == 1990)
df_autos_2000t <- df_autos_clasiftot %>% filter(decada == 2000)

x1960t <- df_autos_1960t$potencia
y1960t <- df_autos_1960t$velmax
x1970t <- df_autos_1970t$potencia
y1970t <- df_autos_1970t$velmax
x1980t <- df_autos_1980t$potencia
y1980t <- df_autos_1980t$velmax
x1990t <- df_autos_1990t$potencia
y1990t <- df_autos_1990t$velmax
x2000t <- df_autos_2000t$potencia
y2000t <- df_autos_2000t$velmax

fit1960t  <- lm(y1960t~log(x1960t))
fit1970t  <- lm(y1970t~log(x1970t))
fit1980t  <- lm(y1980t~log(x1980t))
fit1990t  <- lm(y1990t~log(x1990t))
fit2000t  <- lm(y2000t~log(x2000t))

xx <- seq(50,180, length=10)

Volvemos a dibujar y vemos:

plot(x1960t,y1960t, col = "green", xlab = "", ylab = "", ylim = c(100,200))
points(x1970t,y1970t, col = "blue")
points(x1980t,y1980t, col = "grey")
points(x1990t,y1990t, col = "red")
points(x2000t,y2000t, col = "orange")
lines(xx, predict(fit1960t, data.frame(x1960t=xx)), col="green")
lines(xx, predict(fit1970t, data.frame(x1970t=xx)), col="blue")
lines(xx, predict(fit1980t, data.frame(x1980t=xx)), col="grey")
lines(xx, predict(fit1990t, data.frame(x1990t=xx)), col="red")
lines(xx, predict(fit2000t, data.frame(x2000t=xx)), col="orange")
title(xlab = "Potencia (hp)", ylab = "Velocidad máxima (km/h)")

Si bien el modelo sigue ajustando el gráfico queda “sucio” y no nos aporta información adicional para este trabajo que es observar tendencias. Podemos sacarle los puntos pero perderíamos las referencias visuales para ver si la curva sigue los datos.

Podemos también comprobar el ajuste de cada modelo consultando las características de cada uno:

summary(fit1960)
## 
## Call:
## lm(formula = y1960 ~ log(x1960))
## 
## Residuals:
##       1       2       3       4       5       6 
## -0.7544 -0.8768  4.0941 -2.4623 -1.0689  1.0682 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -3.991     12.416  -0.321 0.763951    
## log(x1960)    32.319      2.745  11.773 0.000298 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.571 on 4 degrees of freedom
## Multiple R-squared:  0.972,  Adjusted R-squared:  0.9649 
## F-statistic: 138.6 on 1 and 4 DF,  p-value: 0.0002978
summary(fit1960t)
## 
## Call:
## lm(formula = y1960t ~ log(x1960t))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.511  -7.400  -0.092   9.544  34.916 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   21.354     29.482   0.724 0.473975    
## log(x1960t)   26.661      6.454   4.131 0.000232 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.62 on 33 degrees of freedom
## Multiple R-squared:  0.3409, Adjusted R-squared:  0.3209 
## F-statistic: 17.07 on 1 and 33 DF,  p-value: 0.0002316

La diferencia entre modelar contra datos preparados vs datos “crudos” puede ser validar la hipótesis o no. Pero también se corre el riesgo de dar por válido algo que no lo es.

Hasta ahora hicimos todo el análisis clasificando los datos según la década. Para algunas cosas esto puede servir pero perdemos visión fina de algunos fenómenos. Estamos metiendo en la misma categoría un auto de 1990 y uno de 1999 por ejemplo. Puede ser que ahí nos estemos perdiendo alguna especificidad. ¿Cómo podemos mostrar esto? Un gráfico común con 50 colores sería complicado de leer y un desafío a la agudeza visual. Lo importante es entender cuando ya no nos sirve un set de herramientas para poder pasar a otras que si sean útiles.

Igual antes un gráfico más de los anteriores pero con el doble de categorías, cada cinco años (corte quinquenal).

df_autos_clasifq <- df_autos_clasiftot %>% mutate(quinq = anio-anio%%5)

df_autos_1960q <- df_autos_clasifq %>% filter(quinq == 1960)
df_autos_1965q <- df_autos_clasifq %>% filter(quinq == 1965)
df_autos_1970q <- df_autos_clasifq %>% filter(quinq == 1970)
df_autos_1975q <- df_autos_clasifq %>% filter(quinq == 1975)
df_autos_1980q <- df_autos_clasifq %>% filter(quinq == 1980)
df_autos_1985q <- df_autos_clasifq %>% filter(quinq == 1985)
df_autos_1990q <- df_autos_clasifq %>% filter(quinq == 1990)
df_autos_1995q <- df_autos_clasifq %>% filter(quinq == 1995)
df_autos_2000q <- df_autos_clasifq %>% filter(quinq == 2000)
df_autos_2005q <- df_autos_clasifq %>% filter(quinq == 2005)

x1960q <- df_autos_1960q$potencia
y1960q <- df_autos_1960q$velmax
x1965q <- df_autos_1965q$potencia
y1965q <- df_autos_1965q$velmax
x1970q <- df_autos_1970q$potencia
y1970q <- df_autos_1970q$velmax
x1975q <- df_autos_1975q$potencia
y1975q <- df_autos_1975q$velmax
x1980q <- df_autos_1980q$potencia
y1980q <- df_autos_1980q$velmax
x1985q <- df_autos_1985q$potencia
y1985q <- df_autos_1985q$velmax
x1990q <- df_autos_1990q$potencia
y1990q <- df_autos_1990q$velmax
x1995q <- df_autos_1995q$potencia
y1995q <- df_autos_1995q$velmax
x2000q <- df_autos_2000q$potencia
y2000q <- df_autos_2000q$velmax
x2005q <- df_autos_2005q$potencia
y2005q <- df_autos_2005q$velmax


fit1960q  <- lm(y1960q~log(x1960q))
fit1965q  <- lm(y1965q~log(x1965q))
fit1970q  <- lm(y1970q~log(x1970q))
fit1975q  <- lm(y1975q~log(x1975q))
fit1980q  <- lm(y1980q~log(x1980q))
fit1985q  <- lm(y1985q~log(x1985q))
fit1990q  <- lm(y1990q~log(x1990q))
fit1995q  <- lm(y1995q~log(x1995q))
fit2000q  <- lm(y2000q~log(x2000q))
fit2005q  <- lm(y2005q~log(x2005q))

xx <- seq(50,180, length=10)

Veamos

plot(x1960q,y1960q, col = "green", xlab = "", ylab = "", ylim = c(100,200))
points(x1965q,y1965q, col = "green")
points(x1970q,y1970q, col = "blue")
points(x1975q,y1975q, col = "blue")
points(x1980q,y1980q, col = "grey")
points(x1985q,y1985q, col = "grey")
points(x1990q,y1990q, col = "red")
points(x1995q,y1995q, col = "red")
points(x2000q,y2000q, col = "orange")
points(x2005q,y2005q, col = "orange")
lines(xx, predict(fit1960q, data.frame(x1960q=xx)), col="green")
lines(xx, predict(fit1965q, data.frame(x1965q=xx)), col="green")
lines(xx, predict(fit1970q, data.frame(x1970q=xx)), col="blue")
lines(xx, predict(fit1975q, data.frame(x1975q=xx)), col="blue")
lines(xx, predict(fit1980q, data.frame(x1980q=xx)), col="grey")
lines(xx, predict(fit1985q, data.frame(x1985q=xx)), col="grey")
lines(xx, predict(fit1990q, data.frame(x1990q=xx)), col="red")
lines(xx, predict(fit1995q, data.frame(x1995q=xx)), col="red")
lines(xx, predict(fit2000q, data.frame(x2000q=xx)), col="orange")
lines(xx, predict(fit2005q, data.frame(x2005q=xx)), col="orange")
title(xlab = "Potencia (hp)", ylab = "Velocidad máxima (km/h)")

Pese a algún corrimiento no hay diferencias significativas dentro de la misma década (líneas del mismo color). Igual es prácticamente ilegible, comparar con el mismo de cinco categorías (o sea cortado por décadas) y se nota ensegiuda.

Vamos a hacer un gráfico de superficie simple. Es la misma matriz de datos de 2x2 pero con la proyección vertical correspondiente. En este caso tendríamos como ejes x e y potencia y año mientras que la elevación (z) sería la velocidad.

Con la librería Akima es muy fácil hacer una interpolación bilinear simple. Primero un poco de preparación para no olvidarse que es cada cosa.

df_autos_m <- df_autos_4 %>% mutate(anio = listaanio - 1900) %>% drop_na() %>% select(-c(listaanio)) %>% filter(listapotencia >= 50 & listapotencia < 150)
df_autos_m$x <- df_autos_m$listapotencia 
df_autos_m$y <- df_autos_m$anio
df_autos_m$z <- df_autos_m$listavel
df_autos_m$listapotencia <- NULL 
df_autos_m$anio <- NULL
df_autos_m$listavel <- NULL

Simplemente hay que indicar que haga la interpolación sobre el rango de datos.

df_matriz <- interp(
    x = df_autos_m$x,
    y = df_autos_m$y,
    z = df_autos_m$z,
    xo = 50:150, 
    yo = 60:110, 
    linear = TRUE,
    duplicate = "strip"
)

Y en el z está la matriz ya armada. En el head salen sólo las primeras seis filas por eso van a ver muchos NA pero funciona.

head(df_matriz$z)
##      [,1] [,2] [,3]  [,4]  [,5]     [,6]     [,7]     [,8]     [,9]
## [1,]   NA   NA  127 125.0 123.0       NA       NA       NA       NA
## [2,]   NA   NA  125 124.4 124.6 126.2069 125.9655 125.7241 125.4828
## [3,]   NA   NA  123 123.8 126.2 129.6552 129.4138 129.1724 128.9310
## [4,]   NA   NA  121 123.2 127.8 133.1034 132.8621 132.6207 131.8163
## [5,]   NA   NA  119 122.6 129.4 136.5517 136.3103 135.2245 133.8571
## [6,]   NA   NA  117 122.0 131.0 140.0000 138.6327 137.2653 135.8980
##         [,10]    [,11]    [,12]    [,13]    [,14]    [,15]    [,16]
## [1,]       NA       NA       NA       NA       NA       NA       NA
## [2,] 125.2414 125.0000 125.8571 126.7143       NA       NA       NA
## [3,] 128.4082 127.0408 127.0000 127.3529 127.7059 128.0588 128.4118
## [4,] 130.4490 129.0816 129.3333 129.0000 129.3529 129.7059 130.0588
## [5,] 132.4898 131.1224 131.6667 131.3333 131.0000 131.3529 131.7059
## [6,] 134.5306 133.1633 133.3889 133.6667 133.3333 133.0000 133.3529
##         [,17]    [,18]    [,19]    [,20]    [,21]    [,22]    [,23]
## [1,]       NA       NA       NA       NA       NA       NA       NA
## [2,]       NA       NA       NA       NA       NA       NA       NA
## [3,] 128.7647 129.1176 129.4706 129.8235 130.4286 133.4000       NA
## [4,] 130.4118 130.7647 131.1176 131.4706 131.8235 132.1765 132.5294
## [5,] 132.0588 132.4118 132.7647 133.1176 133.4706 133.8235 134.1765
## [6,] 133.7059 134.0588 134.4118 134.7647 135.1176 135.4706 135.8235
##         [,24]    [,25]    [,26]    [,27]    [,28]    [,29]  [,30]    [,31]
## [1,]       NA       NA       NA       NA       NA       NA     NA       NA
## [2,]       NA       NA       NA       NA       NA       NA     NA       NA
## [3,]       NA       NA       NA       NA       NA       NA     NA       NA
## [4,] 132.8824 133.2353 133.5882 133.9412 134.2941 134.6471 135.00 142.2000
## [5,] 134.5294 134.8824 135.2353 135.5882 135.9412 137.9767 137.75 136.1250
## [6,] 136.1765 136.5294 136.8824 137.2353 140.9535 145.2326 140.50 136.1042
##         [,32]    [,33]    [,34]    [,35]    [,36] [,37]    [,38]    [,39]
## [1,]       NA       NA       NA       NA       NA    NA       NA       NA
## [2,]       NA       NA       NA       NA       NA    NA       NA       NA
## [3,]       NA       NA       NA       NA       NA    NA       NA       NA
## [4,]       NA       NA       NA       NA       NA    NA       NA       NA
## [5,] 137.2708 138.4167 139.5625 140.7083 141.8542 143.0 145.6667 148.3333
## [6,] 137.2500 138.3958 139.5417 140.6875 141.8333 144.4 146.3333 149.0000
##      [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## [1,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [2,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [3,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [4,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [5,]   151    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [6,]   151    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
##      [,51]
## [1,]    NA
## [2,]    NA
## [3,]    NA
## [4,]    NA
## [5,]    NA
## [6,]    NA

Lo pasamos a una matriz y graficamos

matriz_autos <- df_matriz$z

plot_ly(z = matriz_autos) %>% add_surface()

No está mal pero seguimos en el mismo problema: le agrega mucha complejidad a algo que debería ser simple de ver. Volvamos a la versión dividida por décadas.

df_autos_prem <- df_autos_clasifpot %>% select(decada, promvel, prompot) %>% transmute(
  x = round(prompot/10),
  y = round((decada-1900)/10),
  z = round(promvel)
)

df_matrizm <- interp(
    x = df_autos_prem$x,
    y = df_autos_prem$y,
    z = df_autos_prem$z,
    xo = 5:15, 
    yo = 6:10, 
    linear = TRUE,
    duplicate = "strip"
)

matriz_autos_2 <- df_matrizm$z

plot_ly(z = matriz_autos_2) %>% add_surface()

Ya casi estamos. Faltan ajustar algunos (varios) detalles. Como está ahora tenemos en el eje X la década (contando desde 1960 como 0, 1970 1 y así hasta 2000 como 4), en el Y la potencia (50 = 0 y así hasta 150 = 10) y por último en el Z la velocidad.

Asignamos los nombres que corresponden y los títulos, ubicamos la cámara para que nos quede en el ángulo solicitado y listo.

colnames(matriz_autos_2) <- c("1960", "1970", "1980", "1990", "2000")
rownames(matriz_autos_2) <- c("50", "60", "70", "80", "90", "100", "110", "120", "130", "140", "150")

axx <- list(
  title = 'Año'
)

axy <- list(

  title = 'Potencia (CV)'
)

axz <- list(
  title = 'Velocidad máxima (km/h)'
)

plot_ly(z = matriz_autos_2, x = colnames(matriz_autos_2), y = rownames(matriz_autos_2)) %>% add_surface(opacity = 0.9) %>% layout(scene = list(xaxis=axx,yaxis=axy,zaxis=axz,aspectmode = "manual", aspectratio = list(x=1, y=1, z=0.8), camera = list(eye = list(x = -1, y = -1, z = 1))))