Uma das maneiras mais simples de se aprender qualquer coisa – inclusive R – é ganhar confiança replicando um exemplo, ainda que você não tenha condições, no início, de entender todas as implicações do que está fazendo. Funciona bem com tudo na vida (ou você sabia andar de bicicleta antes de subir nela pela primeira vez?) e, claro, também no R.
Depois de replicar, claro, é legal mudar a base de dados e treinar novamente. Assim, após replicar o exemplo de cross-validation do BlogR, resolvi, como economista, fazer um exercício de, digamos, estática comparativa. Tudo o mais constante, mudei a base de dados. Mas o que vamos fazer? O que é cross-validation, ou melhor, o que é leave-one-out cross-validation? Citando o autor:
Leave-one-out is a type of cross validation whereby the following is done for each observation in the data:
- Run model on all other observations
- Use model to predict value for observation
This means that a model is fitted, and a predicted is made n times where n is the number of observations in your data.
Seu exemplo usa a (já) clássica base de dados do R, mtcars. Eu resolvi usar uma do pacote AER , a base OECDGrowth.
Um ponto importante é que a amostra precisa ser formatada para uso com os comandos do exemplo original. Veremos isso a seguir.
A regressão que vou replicar é uma equação de crescimento que Kleiber e Zeileis estimam na seção 4.4 de seu ótimo livro (obviamente, o Applied Econometrics with R), quando discutem resistant regressions. Trata-se desta equação:
solow_lm <- lm(log(gdp85/gdp60) ~ log(gdp60) + log(invest) + log(popgrowth + .05), data = OECDGrowth)
Vejamos o que fiz.
Primeiro criei as variáveis já em logaritmos para evitar possíveis problemas com o uso de fórmulas dentro das funções utilizadas no exemplo (por via das dúvidas, prefiro me precaver). Eis o código.
library(AER)
data("OECDGrowth")
# vamos criar a base das variaveis usadas na solow_lm
solow_lm <<- lm(log(gdp85/gdp60) ~ log(gdp60) +
+ log(invest) + log(popgrowth + .05), data = OECDGrowth)
gdp85<-OECDGrowth$gdp85
gdp60<-OECDGrowth$gdp60
invest<-OECDGrowth$invest
popgrowth<-OECDGrowth$popgrowth
lgdp85_60<-log(gdp85/gdp60)
lgdp60<-log(gdp60)
linves<-log(invest)
lpopgrowthplus005<-log(popgrowth+0.05)
data1<-data.frame(lgdp85_60,lgdp60,linvest,lpopgrowthplus005)
Pronto para a replicação? É bom ler o código original com atenção. Ao fazer isso, descobre-se facilmente que é preciso alterar o nome da base de dados original – e o da variável dependente – para não ter surpresas desagradáveis. O trecho do código abaixo pode ser comparado com o original citado lá no início deste pequeno texto.
library(pipelearner)
pl <- pipelearner(data1, lm, lgdp85_60 ~ .)
pl <- learn_cvpairs(pl, k = nrow(data1))
pl <- learn(pl)
# ou
pl <- pipelearner(data1, lm, lgdp85_60 ~ .) %>%
learn_cvpairs(k = nrow(data1)) %>%
learn()
O que está acontecendo? Não é esta a minha preocupação, neste estágio do exercício (como disse no início do texto). Mas você pode descobrir um pouco mais aqui. Basicamente, estamos usando um algoritmo de machine learning na regressão de crescimento original.
Os passos seguintes não necessitam alterações. Para não ficar tão igual, vou traduzir os elementos textuais do gráfico (ok, eu sei que é preguiça).
library(tidyverse)
# Extrai os valores observados e previstos da variavel dependente para cada observação
pl <- pl %>%
mutate(true = map2_dbl(test, target, ~as.data.frame(.x)[[.y]]),
predicted = map2_dbl(fit, test, predict))
# Resumo dos resultados
results <- pl %>%
summarise(
sse = sum((predicted - true)^2),
sst = sum(true^2)
) %>%
mutate(r_squared = 1 - sse / sst)
results
pl %>%
ggplot(aes(true, predicted)) +
geom_point(size = 2) +
geom_abline(intercept = 0, slope = 1, linetype = 2) +
theme_minimal() +
labs(x = "Valor observado", y = "Valor previsto") +
ggtitle("Valores Observados contra previstos baseados\n em leave-one-one cross validation")
Agora, vejamos o gráfico.
Comparando com o gráfico do autor para seu exercício em mtcars confesso que não consigo dizer se há alguma superioridade no ajuste. Isto é bastante razoável, dado que o R² obtido por ele (0.95) não é tão distante do que o obtido aqui (0.94).
Ok, concluímos o nosso objetivo de refazer os passos do exercício original com outra base de dados. Vale ressaltar que o pipelearner tem funções mais interessantes e talvez você queira estudar outros exemplos criados pelo autor do pacote (aliás, caso você queira realmente entender o que foi feito aqui, é bom dar uma olhada lá). Aliás, após o sucesso de um exercício como este, a gente ganha confiança para, justamente, seguir em busca de maior compreensão, não é?
Até a próxima.
