Introdução

Estamos interessados em explicar o comportamento da bilheteria de filmes exibidos no cinema, para isso devemos propor um modelo preditivo de regressão/machine learning. Temos um conjunto de dados com 2.400 filmes e variáveis como: Orçamento, gêneros, título, popularidade, compania de produção, país de produção, data de estréia, duração, idiomas, atores e bilheteria. Algumas variáveis estão com muitas informações em um campo, portanto iremos filtrar algumas variáveis. Para nosso estudo decidimos adotar as variáveis: Orçamento, popularidade, duração e país para explicar a variável bilheteria (revenue).

library(rsample)
library(gbm)
library(xgboost)
library(caret)
library(h2o)
library(pdp)
library(ggplot2)
library(lime)
library(tidyverse)
library(readr)
set.seed(519)

train <- train %>%
  mutate(
    Pais = str_extract(
      production_countries, pattern = "\\[\\{'iso_3166_1': '[\\w\\d ]*',"),
    Pais = str_sub(Pais, start = nchar(Pais) - 3, end = nchar(Pais) - 2)
  ) %>% select("budget", "popularity", "runtime", "revenue", "Pais")


true <- true %>%
  mutate(
    Pais = str_extract(
      production_countries, pattern = "\\[\\{'iso_3166_1': '[\\w\\d ]*',"),
    Pais = str_sub(Pais, start = nchar(Pais) - 3, end = nchar(Pais) - 2)
  ) %>% select("budget", "popularity", "runtime", "revenue", "Pais")

train$Pais <- factor(train$Pais)
true$Pais <- factor(true$Pais)

GBM: Gradient Boosting Machine

Para a construção do nosso modelo utilizaremos o pacote GBM (Gradient Boosting Machine). Começaremos ajustando um modelo gbm com a configuração padrão, tendo assim um resultado inicial.

gbm.fit <- gbm(
  formula = revenue ~ .,
  distribution = "gaussian",
  data = train,
  n.trees = 5000,
  interaction.depth = 1,
  shrinkage = 0.001,
  cv.folds = 5,
  n.cores = NULL,
  verbose = FALSE
)

#resultados
print(gbm.fit)
## gbm(formula = revenue ~ ., distribution = "gaussian", data = train, 
##     n.trees = 5000, interaction.depth = 1, shrinkage = 0.001, 
##     cv.folds = 5, verbose = FALSE, n.cores = NULL)
## A gradient boosted model with gaussian loss function.
## 5000 iterations were performed.
## The best cross-validation iteration was 5000.
## There were 4 predictors of which 4 had non-zero influence.
#pega o RMSE da melhor iteração
sqrt(min(gbm.fit$cv.error))
## [1] 83393188
gbm.perf(gbm.fit, method = "cv")

## [1] 5000

Agora vamos “afinar” o modelo (tuning), isto é, procurar os parâmetros que nos fornece o melhor modelo. Poderíamos alterar os parâmetros manualmente de forma a procurar os ideais, mas optaremos por fazer de forma automática. Faremos uma grid que irá gerar todas as combinações possíveis de alguns valores para os parâmetros. Como podemos ver, temos 16 combinações de parâmetros.

#grid 1
hyper_grid <- expand.grid(
  shrinkage = c(.001, .1, .3),
  interaction.depth = c(1, 3, 5),
  n.minobsinnode = c(5, 10, 15),
  bag.fraction = c( .8),
  optimal_trees = 0,               
  min_RMSE = 0                    
)
nrow(hyper_grid)
## [1] 27

Agora iremos computar os 16 possíveis modelos de gbm com as combinações de parâmetros geradas anteriormente.

for(i in 1:nrow(hyper_grid)) {
  
  gbm.tune <- gbm(
    formula = revenue ~ .,
    distribution = "gaussian",
    data = train,
    n.trees = 5000,
    interaction.depth = hyper_grid$interaction.depth[i],
    shrinkage = hyper_grid$shrinkage[i],
    n.minobsinnode = hyper_grid$n.minobsinnode[i],
    bag.fraction = hyper_grid$bag.fraction[i],
    train.fraction = .8,
    n.cores = NULL, # will use all cores by default
    verbose = FALSE
  )

  hyper_grid$optimal_trees[i] <- which.min(gbm.tune$valid.error)
  hyper_grid$min_RMSE[i] <- sqrt(min(gbm.tune$valid.error))
}

hyper_grid %>%
  dplyr::arrange(min_RMSE) %>%
  head(10)
##    shrinkage interaction.depth n.minobsinnode bag.fraction optimal_trees
## 1      0.100                 5              5          0.8            27
## 2      0.300                 5              5          0.8             9
## 3      0.100                 3              5          0.8            55
## 4      0.100                 5             10          0.8            34
## 5      0.300                 3              5          0.8            14
## 6      0.100                 3             10          0.8            42
## 7      0.001                 5             10          0.8          3026
## 8      0.001                 5             15          0.8          2720
## 9      0.100                 5             15          0.8            30
## 10     0.001                 5              5          0.8          3421
##    min_RMSE
## 1  83124746
## 2  83614165
## 3  83788505
## 4  83924257
## 5  83956045
## 6  84162749
## 7  84317889
## 8  84355639
## 9  84409103
## 10 84659946

Portanto o menor RMSE possível com essas combinações de parâmetros é de 83124746 que já é menos que o do modelo inicial, portanto já temos um modelo melhor do que o inicial, mas vamos tentar encontrar um ainda melhor. Portanto procuraremos por outras combinações, de forma que elas estejam próximas dos valores do parâmetros que tiveram menor RMSE na grid anterior. O shrinkage com melhores resultados foi de 0,1 e 0,3, portanto buscaremos por 0,1 ; 0,2 e 0,3. Faremos escolhas usando a mesma lógica para os outros parâmetros também.

#grid 2
hyper_grid <- expand.grid(
  shrinkage = c( .1, .2, .3),
  interaction.depth = c(5, 7),
  n.minobsinnode = c(5, 10, 15),
  bag.fraction = c(.8),
  optimal_trees = 0,
  min_RMSE = 0
)
nrow(hyper_grid)
## [1] 18
for(i in 1:nrow(hyper_grid)) {

  gbm.tune <- gbm(
    formula = revenue ~ .,
    distribution = "gaussian",
    data = train,
    n.trees = 5000,
    interaction.depth = hyper_grid$interaction.depth[i],
    shrinkage = hyper_grid$shrinkage[i],
    n.minobsinnode = hyper_grid$n.minobsinnode[i],
    bag.fraction = hyper_grid$bag.fraction[i],
    train.fraction = .8,
    n.cores = NULL, # will use all cores by default
    verbose = FALSE
  )

  hyper_grid$optimal_trees[i] <- which.min(gbm.tune$valid.error)
  hyper_grid$min_RMSE[i] <- sqrt(min(gbm.tune$valid.error))
}

hyper_grid %>%
  dplyr::arrange(min_RMSE) %>%
  head(10)
##    shrinkage interaction.depth n.minobsinnode bag.fraction optimal_trees
## 1        0.2                 7             15          0.8            15
## 2        0.2                 5              5          0.8            12
## 3        0.2                 7             10          0.8            22
## 4        0.1                 5             10          0.8            34
## 5        0.2                 5             10          0.8            12
## 6        0.3                 7             10          0.8             6
## 7        0.1                 7             15          0.8            37
## 8        0.3                 5             15          0.8            16
## 9        0.1                 7             10          0.8            31
## 10       0.3                 7              5          0.8            14
##    min_RMSE
## 1  83296229
## 2  83312145
## 3  83754402
## 4  83981655
## 5  84023559
## 6  84513291
## 7  84656921
## 8  84752343
## 9  84783004
## 10 84926006

Agora o menor RMSE é de 83296229, que é maior que o encontrado no grid anterior. Portanto ficaremos com o anterior, sendo então o nosso modelo final.

gbm.fit <- gbm(
  formula = revenue ~ .,
  distribution = "gaussian",
  data = train,
  n.trees = 100,
  interaction.depth = 5,
  shrinkage = 0.1,
  cv.folds = 5,
  n.cores = NULL,
  verbose = FALSE
)

#resultados
print(gbm.fit)
## gbm(formula = revenue ~ ., distribution = "gaussian", data = train, 
##     n.trees = 100, interaction.depth = 5, shrinkage = 0.1, cv.folds = 5, 
##     verbose = FALSE, n.cores = NULL)
## A gradient boosted model with gaussian loss function.
## 100 iterations were performed.
## The best cross-validation iteration was 28.
## There were 4 predictors of which 4 had non-zero influence.
#pega o RMSE da melhor iteração
sqrt(min(gbm.fit$cv.error))
## [1] 81763864
gbm.perf(gbm.fit, method = "cv")

## [1] 28

Sendo esse o nosso modelo final. Podemos verificar a seguir quais variáveis explicam melhor/influenciam mais na bilheteria dos filmes. Como esperado o orçamento de um filme é o que mais influencia em sua bilheteria, a popularidade dos filmes também influencia razoávelmente bem. Já a duração e país do filme não tem muita influência na bilheteria.

par(mar = c(5, 8, 1, 1))
summary(
  gbm.fit,
  cBars = 10,
  method = relative.influence, # also can use permutation.test.gbm
  las = 2
)

##                   var   rel.inf
## budget         budget 62.467177
## popularity popularity 26.226661
## runtime       runtime  6.815979
## Pais             Pais  4.490183