map2 (funções que recebem dois argumentos)

Em diversos cenários precisamos trabalhar com funções com mais de um argumento. Considere um exemplo em que queremos calcular a área de diferentes retângulos a partir de um vetor com as bases e outro com as alturas. Poderíamos utilizar a função map2 da seguinte forma:

base <- c(1, 2, 3)
altura <- c(4, 5, 6)

map2(base, altura, prod)

## [[1]]
## [1] 4
## 
## [[2]]
## [1] 10
## 
## [[3]]
## [1] 18

A partir de uma estrutura de banco de dados, considere que queremos obter a soma das variáveis num1 e num2 para cada elemento. Aqui consideraremos a função map2 e map2_int.

tibble(num1 = 1:5, 
       num2 = -15:-11) %>% 
  mutate(soma = map2(num1, num2, sum))

## # A tibble: 5 x 3
##    num1  num2 soma     
##   <int> <int> <list>   
## 1     1   -15 <int [1]>
## 2     2   -14 <int [1]>
## 3     3   -13 <int [1]>
## 4     4   -12 <int [1]>
## 5     5   -11 <int [1]>

tibble(num1 = 1:5, 
       num2 = -15:-11) %>% 
  mutate(soma = map2_int(num1, num2, sum))

## # A tibble: 5 x 3
##    num1  num2  soma
##   <int> <int> <int>
## 1     1   -15   -14
## 2     2   -14   -12
## 3     3   -13   -10
## 4     4   -12    -8
## 5     5   -11    -6

Note que podemos fazer isso de forma diferente, mas o objetivo aqui é facilitar a compreensão e aplicação das funções da família map.

Quando a função precisa receber 3 ou mais argumentos, podemos utilizar a função pmap.

pmap (funções que recebem 3 ou mais argumentos)

Considere o cenário em que queremos obter, para cada elemento do banco de dados, a soma dos valores de três colunas (num1, num2 e num3). Podemos utilizar as funções pmap e pmap_dbl.

tibble(num1 = 1:5, 
       num2 = -15:-11, 
       num3 = seq(10, 50, 10)) %>% 
  mutate(soma = pmap(list(num1, num2, num3), sum))

## # A tibble: 5 x 4
##    num1  num2  num3 soma     
##   <int> <int> <dbl> <list>   
## 1     1   -15    10 <dbl [1]>
## 2     2   -14    20 <dbl [1]>
## 3     3   -13    30 <dbl [1]>
## 4     4   -12    40 <dbl [1]>
## 5     5   -11    50 <dbl [1]>

tibble(num1 = 1:5, 
       num2 = -15:-11, 
       num3 = seq(10, 50, 10)) %>% 
  mutate(soma = pmap_dbl(list(num1, num2, num3), sum))

## # A tibble: 5 x 4
##    num1  num2  num3  soma
##   <int> <int> <dbl> <dbl>
## 1     1   -15    10    -4
## 2     2   -14    20     8
## 3     3   -13    30    20
## 4     4   -12    40    32
## 5     5   -11    50    44

Em muitas situações, criaremos uma função para um único processamento. Para evitar a criação de um objeto desnecessário, utilizaremos a estrutura de fórmula que será convertida em uma função diretamente dentro do map. Nesse caso, temos três formas diferentes de nos referir aos argumentos:

tibble(numero = c(1, 4, 9, 16, 25)) %>% 
  mutate(raiz = map_dbl(numero, ~ .x + 2))

## # A tibble: 5 x 2
##   numero  raiz
##    <dbl> <dbl>
## 1      1     3
## 2      4     6
## 3      9    11
## 4     16    18
## 5     25    27

tibble(numero = c(1, 4, 9, 16, 25), 
       valor = c(10, 3, 4, 5, 8)) %>% 
  mutate(raiz = map2_dbl(numero, valor,  ~ .x + .y))

## # A tibble: 5 x 3
##   numero valor  raiz
##    <dbl> <dbl> <dbl>
## 1      1    10    11
## 2      4     3     7
## 3      9     4    13
## 4     16     5    21
## 5     25     8    33

tibble(numero = c(1, 4, 9, 16, 25), 
       valor = c(10, 3, 4, 5, 8), 
       fator = c(.1, .5, .9, .7, .2)) %>% 
  mutate(raiz = pmap_dbl(list(numero, valor, fator),  ~ ..3 * (..1 + ..2)))

## # A tibble: 5 x 4
##   numero valor fator  raiz
##    <dbl> <dbl> <dbl> <dbl>
## 1      1    10   0.1   1.1
## 2      4     3   0.5   3.5
## 3      9     4   0.9  11.7
## 4     16     5   0.7  14.7
## 5     25     8   0.2   6.6

Ajuste de hiperparâmetro por validação cruzada

Nessa aplicação, como base, utilizaremos o conjunto de dados Credit do pacote ISLR. O objetivo desse banco de dados é predizer os valores da variável Balance.

dados <- ISLR::Credit 

head(dados)

Suponha que será utilizado um modelo de Floresta Aleatória e deseja-se definir o hiperparâmetro mtry (define o número de variáveis que serão sorteadas para cada divisão dos nós). Queremos testar os valores 2, 4, 8 e 10 para essa quantidade e, consequentemente, estimar qual desses valores leva a um menor erro de previsão. Para isso, utilizaremos a validação cruzada. A função vold_cv nos auxiliará nesse processo (veja o post Descartes e rsample para mais detalhes sobre a função vfold_cv).

set.seed(15)

vfold_cv(dados, v = 10) # v indica o número de lotes

## #  10-fold cross-validation 
## # A tibble: 10 x 2
##    splits           id    
##    <list>           <chr> 
##  1 <split [360/40]> Fold01
##  2 <split [360/40]> Fold02
##  3 <split [360/40]> Fold03
##  4 <split [360/40]> Fold04
##  5 <split [360/40]> Fold05
##  6 <split [360/40]> Fold06
##  7 <split [360/40]> Fold07
##  8 <split [360/40]> Fold08
##  9 <split [360/40]> Fold09
## 10 <split [360/40]> Fold10

Para cada uma das dez linhas são definidos conjuntos de treinamento (com 360 observações) e conjuntos de teste (40 observações). O lote é identificado pela coluna id. Faremos a combinação desses lotes com os valores experimentais de mtry.

vfold_cv(dados, v = 10) %>% 
  crossing(mtry = c(2, 4, 8, 10))

## # A tibble: 40 x 3
##    splits           id      mtry
##    <list>           <chr>  <dbl>
##  1 <split [360/40]> Fold01     2
##  2 <split [360/40]> Fold01     4
##  3 <split [360/40]> Fold01     8
##  4 <split [360/40]> Fold01    10
##  5 <split [360/40]> Fold02     2
##  6 <split [360/40]> Fold02     4
##  7 <split [360/40]> Fold02     8
##  8 <split [360/40]> Fold02    10
##  9 <split [360/40]> Fold03     2
## 10 <split [360/40]> Fold03     4
## # ... with 30 more rows

Assim, para cada linha do tibble acima, ajustaremos a Floresta Aleatória com o conjunto de treino do split considerando o mtry da respectiva linha e avaliaremos o erro no conjunto de teste. Isso será feito com a aplicação de uma função para cada linha do tibble com map_dbl.

Primeiro, precisaremos definir a função para calcular o erro quadrático médio nessa situação.

ajuste <- function(splits, variaveis_mtry) {
  
  treino <- training(splits)
  teste <- testing(splits)
  
  rf <- ranger(Balance ~ . - ID, mtry = variaveis_mtry, data = treino)
  
  eqm <- mean((teste$Balance - predict(rf, teste)$predictions)^2)
  
  return(eqm)
  
}

Pronto! Agora, para obter o erro estimado em cada cenário, aplicaremos a função ajustecom auxílio da função map (execute o comando a seguir linha por linha e verifique as saídas).

set.seed(321)

vfold_cv(dados, v = 10) %>% 
  crossing(mtry = c(2, 4, 8, 10)) %>% 
  mutate(eqm = map2_dbl(splits, mtry, ajuste)) %>% 
  group_by(mtry) %>% 
  summarise(eqm_medio = mean(eqm), 
            desv_pad = sd(eqm)) %>% 
  arrange(eqm_medio)

## # A tibble: 4 x 3
##    mtry eqm_medio desv_pad
##   <dbl>     <dbl>    <dbl>
## 1    10     9950.    3959.
## 2     8    10318.    4039.
## 3     4    14760.    5630.
## 4     2    27114.    8353.

Portanto, o valor de mtry que levou a menor estimativa pontual do erro quadrático médio foi de 10. Dessa forma, utilizaremos 10 variáveis em cada split das árvores da Floresta Aleatória.

E se tivéssemos como objetivo avaliar o ajuste do hiperparâmetro mtry e min.node.size? Poderíamos utilizar a função pmap_dlb da seguinte forma:

ajustes2 <- function(splits, variaveis_mtry, node_size) {
  
  treino <- training(splits)
  teste <- testing(splits)
  
  rf <- ranger(Balance ~ . - ID, mtry = variaveis_mtry, min.node.size = node_size, data = treino)
  
  eqm <- mean((teste$Balance - predict(rf, teste)$predictions)^2)
  
  return(eqm)
  
}

set.seed(321)

vfold_cv(dados, v = 10) %>% 
  crossing(mtry = c(2, 4, 8, 10), 
           node_size = c(5, 10, 15)) %>% 
  mutate(eqm = pmap_dbl(list(splits, mtry, node_size), ajustes2)) %>% 
  group_by(mtry, node_size) %>% 
  summarise(eqm_medio = mean(eqm), 
            desv_pad = sd(eqm)) %>% 
  arrange(eqm_medio)

## # A tibble: 12 x 4
## # Groups:   mtry [4]
##     mtry node_size eqm_medio desv_pad
##    <dbl>     <dbl>     <dbl>    <dbl>
##  1    10         5     9876.    3804.
##  2     8         5    10158.    3967.
##  3    10        10    10490.    3926.
##  4     8        10    11347.    4598.
##  5    10        15    11838.    4680.
##  6     8        15    12312.    5127.
##  7     4         5    14523.    5623.
##  8     4        10    15948.    6101.
##  9     4        15    17730.    6611.
## 10     2         5    26583.    8001.
## 11     2        10    29425.    9239.
## 12     2        15    32271.    9782.

Logo, nesse cenário, o ideal seria utilizar mtry = 10 e node_size = 5.

Note como podemos facilitar diversas tarefas apenas com funções do pacote rsample e com a definição de funções simples.

A seguir veremos algumas aplicações para entender um pouco mais sobre a influência de certas quantidades em conceitos estatísticos.

Intervalo de Confiança

Nesta simulação trabalharemos com intervalos de confiança para a média com dados de uma distribuição normal com variância conhecida. Lembre-se que o intervalo de confiança para esse caso é dado por:

\[ IC(\mu, \gamma) = \bigg[\overline{X} \pm z_{(1 + \gamma)/2} \frac{\sigma}{\sqrt{n}}\bigg], \]

em que \(\sigma\) é o desvio padrão, \(n\) é o tamanho da amostra e \(\gamma\) é o nível de confiança.

Inicialmente consideraremos uma amostra (n_amostra) de tamanho 100 de uma população com média (media_pop) igual a 50, desvio padrão (dp_pop) igual a 3 e um nível de confiança (confianca) de 95%. Vamos criar a função simulacao para retornar a média amostral e os limites inferior e superior do intervalo.

media_pop <- 50
dp_pop <- 3
n_amostra <- 100
confianca <- .95


simulacao <- function(n_amostra, media_pop, dp_pop, confianca = .95) {
  
  x <- rnorm(n_amostra, media_pop, dp_pop)
  
  tibble(media_amostral = mean(x), 
         lim_inf = mean(x) - qnorm((1 + confianca)/2) * dp_pop / sqrt(n_amostra),
         lim_sup = mean(x) + qnorm((1 + confianca)/2) * dp_pop / sqrt(n_amostra))
  
}

Agora, aplicaremos a função simulacao 100 vezes. Veja a saída.

tibble(b = 1:10^2) %>% 
  mutate(resultado = map(b, ~simulacao(n_amostra, media_pop, dp_pop, confianca))) %>% 
  unnest(cols = resultado) 

## # A tibble: 100 x 4
##        b media_amostral lim_inf lim_sup
##    <int>          <dbl>   <dbl>   <dbl>
##  1     1           50.4    49.8    51.0
##  2     2           50.2    49.7    50.8
##  3     3           50.2    49.6    50.8
##  4     4           50.1    49.5    50.7
##  5     5           50.1    49.5    50.6
##  6     6           49.9    49.3    50.5
##  7     7           49.5    48.9    50.0
##  8     8           49.6    49.0    50.2
##  9     9           50.5    49.9    51.0
## 10    10           50.1    49.5    50.7
## # ... with 90 more rows

Identificaremos com uma cor vermelha os intervalos que não contenham a média populacional identificada pela linha azul tracejada. Qual a porcentagem de intervalos contendo a média populacional você esperaria?

tibble(b = 1:10^2) %>% 
  mutate(resultado = map(b, ~simulacao(n_amostra, media_pop, dp_pop, confianca))) %>% 
  unnest(cols = resultado) %>% 
  mutate(contem = case_when(media_pop >= lim_inf & media_pop <= lim_sup ~ "contém", 
                            TRUE ~ "não contém")) %>% 
  ggplot(aes(b, media_amostral, color = contem)) + 
    geom_hline(yintercept = media_pop, color = "blue", lty = 2) +
    geom_errorbar(aes(ymin = lim_inf, ymax = lim_sup)) +
    labs(x = "simulação", y = NULL, color = NULL) +
    scale_color_manual(values = c("contém" = "black", "não contém" = "red")) +
    theme_bw() + 
    theme(legend.position = "none")

Intervalo de confiança - variando o nível de confiança

O que aconteceria com os intervalos se considerássemos diferentes níveis de confiança? Um nível de confiança maior levaria a um intervalo de amplitude maior? Para verificar essas questões, faremos uma simulação variando essa quantidade com os conceitos apresentados.

crossing(b = 1:10^2, 
         confianca = c(.25, .50, .75, .99)) %>% 
  mutate(resultado = map(confianca, ~simulacao(n_amostra, media_pop, dp_pop, .x))) %>% 
  unnest(cols = resultado) %>% 
  mutate(contem = case_when(media_pop >= lim_inf & media_pop <= lim_sup ~ "contém", 
                            TRUE ~ "não contém")) %>% 
  ggplot(aes(b, media_amostral, color = contem)) + 
    geom_hline(yintercept = media_pop, color = "blue", lty = 2) +
    geom_errorbar(aes(ymin = lim_inf, ymax = lim_sup)) +
    labs(x = "simulação", y = NULL, color = NULL) +
    facet_wrap(~ confianca, ncol = 2) +
    scale_color_manual(values = c("contém" = "black", "não contém" = "red")) +
    theme_bw() + 
    theme(legend.position = "none")

Note que baixos valores do nível de confiança estão associados a intervalos menores e um número reduzido de intervalos contendo a média populacional. Veja como os intervalos com 99% de confiança apresentam uma grande amplitude.

Intervalo de confiança - variando o tamanho amostral

E o que aconteceria com a amplitude dos intervalos ao variar o tamanho amostral fixando-se o nível de confiança?

crossing(b = 1:10^2, 
         n = c(10, 100, 500, 1000)) %>% 
  mutate(resultado = map(n, ~simulacao(.x, media_pop, dp_pop, confianca))) %>% 
  unnest(cols = resultado) %>% 
  mutate(contem = case_when(media_pop >= lim_inf & media_pop <= lim_sup ~ "contém", 
                            TRUE ~ "não contém")) %>% 
  ggplot(aes(b, media_amostral, color = contem)) + 
    geom_hline(yintercept = media_pop, color = "blue", lty = 2) +
    geom_errorbar(aes(ymin = lim_inf, ymax = lim_sup)) +
    labs(x = "simulação", y = NULL, color = NULL) +
    facet_wrap(~ n, ncol = 2) +
    scale_color_manual(values = c("contém" = "black", "não contém" = "red")) +
    theme_bw() + 
    theme(legend.position = "none")

Intervalo de confiança - variando o tamanho amostral e o nível de confiança

E como se relacionam as variações do tamanho amostral e nível de confiança? A seguir apresentamos uma simulação testando diferentes combinações.

crossing(b = 1:10^2, 
         n = c(100, 500, 1000), 
         confianca = c(.33, .66, .99)) %>% 
  mutate(resultado = map2(n, confianca, ~simulacao(.x, media_pop, dp_pop, .y))) %>% 
  unnest(cols = resultado) %>% 
  mutate(contem = case_when(media_pop >= lim_inf & media_pop <= lim_sup ~ "contém", 
                            TRUE ~ "não contém")) %>% 
  ggplot(aes(b, media_amostral, color = contem)) + 
    geom_hline(yintercept = media_pop, color = "blue", lty = 2) +
    geom_errorbar(aes(ymin = lim_inf, ymax = lim_sup)) +
    labs(x = "simulação", y = NULL, color = NULL) +
    facet_wrap(n ~ confianca, ncol = 3) +
    scale_color_manual(values = c("contém" = "black", "não contém" = "red")) +
    theme_bw() + 
    theme(legend.position = "none")

Conclusão

Nesse post aprendemos a utilizar as funções da família map com algum detalhe. Verificamos como fazer o ajuste de hiperparâmetros de forma muito rápida com auxílio do pacote rsample e como realizar diferentes simulações com poucas linhas de código. Esse último contexto de aplicação é muito útil para quem está começando a estudar estatística.

Caso tenha alguma crítica, sugestão ou comentário, me envie uma mensagem!