Introdução

O objetivo desse estudo é construir um modelo de regressão linear para investigar quais fatores podem influenciar no dinheiro arrecadado com a bilheteria de filmes de cinema. Os dados foram obtidos através de um conjunto de dados disponibilizado no Kaggle, que contém informações como: número de críticas no IMDB, duração do filme, quantidade de curtidas na página do Facebook do diretor, quantidade de curtidas na página do Facebook dos atores, dinheiro arrecadado com a bilheteria, número de usuários que avaliaram no IMDB, censura, nota no IMDB, orçamento do filme. Espera-se construir um modelo para explicar a relação linear entre a bilheteria e essas outras informações.

Análise inicial

Para haver consistência nos dados decidiu-se utilizar apenas os filmes exibidos nos cinemas dos Estados Unidos no ano de 2014, de forma a garantir que variáveis como bilheteria e orçamento estariam em dólares. A amostra obtida tem 120 observações.

O conjunto de dados contava inicialmente com 28 variáveis, portanto algumas são dispensáveis para nossa análise, como: filme colorido ou preto e branco, título do filme, proporção da tela. Após eliminar variáveis como essa sobraram 13 além da que estamos interessados (bilheteria). Em seguida utilizou-se o método de seleção automática stepwise com o critério de informação de Akaike (AIC), o modelo obtido foi:

## 
## Call:
## lm(formula = gross ~ director_facebook_likes + num_voted_users + 
##     budget + movie_facebook_likes, data = dados)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -122097320  -24693444   -7502718   15299542  216953657 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              1.323e+07  6.149e+06   2.151   0.0336 *  
## director_facebook_likes  2.906e+03  1.294e+03   2.246   0.0266 *  
## num_voted_users          2.748e+02  6.143e+01   4.473 1.82e-05 ***
## budget                   7.325e-01  9.890e-02   7.407 2.34e-11 ***
## movie_facebook_likes    -5.269e+02  2.059e+02  -2.559   0.0118 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 45640000 on 115 degrees of freedom
## Multiple R-squared:  0.641,  Adjusted R-squared:  0.6285 
## F-statistic: 51.33 on 4 and 115 DF,  p-value: < 2.2e-16

Em seguida apresenta-se um diagrama com os gráficos de dispersão e coeficientes de correlação para os pares de variáveis que serão utilizadas.

Modelos iniciais propostos

Através do resultado obtido pela utilização do stepAIC anteriormente propõe-se o seguinte modelo. Supondo que \(\epsilon_{i} \sim N(0,\sigma^{2})\) e independentes: \[y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \beta_4 x_{4i} + \epsilon_{i}\] onde \(i = 1, 2, ..., 120\), \(y\) = Bilheteria em dólares, \(x_1\) = Orçamento em dólares, \(x_2\) = Quantidade de avaliações dos usuários, \(x_3\) = Número de curtidas na página do Facebook do diretor, \(x_4\) = Número de curtidas na página do Facebook do filme.

Para garantir a validade do modelo iremos verificar as suposições feitas. Para testar a suposição de normalidade dos erros sugere-se o teste de Bera-Jarque. Considerando um nível de significância \(\alpha = 0,05\), como o p-valor obtido no teste de Bera-Jarque foi muito pequeno rejeita-se a hipótese nula de normalidade dos erros. E para testar a suposição de homoscedasticidade dos erros usaremos o teste de Koenker, visto que o teste de Goldfed-Quandt é muito dependente da normalidade e o teste de Breusch-Pagan costuma ser menos poderoso que o teste de Koenker sob não normalidade. Para o teste de Koenker o p-valor também foi pequeno, portanto rejeita-se a hipótese nula de homoscedasticidade.

## 
##  Jarque Bera Test
## 
## data:  res_AIC
## X-squared = 182.72, df = 2, p-value < 2.2e-16
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste_AIC
## BP = 16.184, df = 4, p-value = 0.002782

Como a bilheteria é sempre positiva utilizaremos a transformação de Box-Cox para reduzir desvios das suposições de normalidade e homoscedasticidade dos erros feitas no modelo, sacrificando interpretabilidade para obter um modelo mais adequado. Primeiro estima-se o \(\lambda\) através da maximização da verossimilhança perfilada.

## [1] 0.4646465

Em seguida realizamos a transformação de Box-Cox com o \(\lambda \approx 0,46\). Como \(\lambda \neq 0\): \[y'(\lambda) = \frac{y^{\lambda} -1}{\lambda}\]

Portanto, ajustando agora um modelo com a variável transformada \(y'\) o teste de Bera-Jarque apresenta um p-valor alto, indicando a não rejeição da hipótese nula de normalidade. O teste de Koenker também apresenta um p-valor alto, indicando a não rejeição da hipótese de homoscedasticidade. Portanto a transformação foi efetiva para reduzir desvios das suposições de normalidade e homoscedasticidade. Porém, nesse novo modelo duas variáveis deixaram de ser significativas. Removendo as variáveis \(x_3\) e \(x_4\) obtemos o modelo final proposto, exibido na próxima seção.

## 
## Call:
## lm(formula = gross_t ~ director_facebook_likes + num_voted_users + 
##     budget + movie_facebook_likes, data = dados_t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7782.9 -1215.3  -124.7  1638.9  8406.0 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              5.202e+03  3.765e+02  13.817  < 2e-16 ***
## director_facebook_likes  1.337e-01  7.923e-02   1.687  0.09432 .  
## num_voted_users          1.206e-02  3.761e-03   3.205  0.00175 ** 
## budget                   4.350e-05  6.056e-06   7.183 7.28e-11 ***
## movie_facebook_likes    -1.911e-02  1.261e-02  -1.516  0.13229    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2795 on 115 degrees of freedom
## Multiple R-squared:  0.5896, Adjusted R-squared:  0.5753 
## F-statistic:  41.3 on 4 and 115 DF,  p-value: < 2.2e-16
## 
##  Jarque Bera Test
## 
## data:  res_t
## X-squared = 1.6145, df = 2, p-value = 0.4461
## 
##  studentized Breusch-Pagan test
## 
## data:  ajuste_t
## BP = 6.4936, df = 4, p-value = 0.1652

Modelo final proposto

Supondo que \(\epsilon_{i} \sim N(0,\sigma^{2})\) e que os erros são independentes: \[y'_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \epsilon_{i}\] onde \(i = 1, 2, ..., 120\), \(y'\) = \(\frac{y^{0,46} -1}{0,46}\), \(x_1\) = Orçamento em dólares, \(x_2\) = Quantidade de avaliações dos usuários.

## 
## Call:
## lm(formula = gross_t ~ budget + num_voted_users, data = dados_t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7914.5 -1523.4  -120.6  1693.7  9578.4 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     5.139e+03  3.759e+02  13.672  < 2e-16 ***
## budget          4.412e-05  5.923e-06   7.449 1.76e-11 ***
## num_voted_users 8.525e-03  2.047e-03   4.165 5.98e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2816 on 117 degrees of freedom
## Multiple R-squared:  0.5761, Adjusted R-squared:  0.5688 
## F-statistic:  79.5 on 2 and 117 DF,  p-value: < 2.2e-16

Ajustando o modelo obtemos:

\[\hat{y'_i} = 5,13.10^{3} + 4,41.10^{-5} x_{1i} + 8,52.10^{-3} x_{2i}\] onde \(i = 1, 2, ..., 120\), \(y'\) = \(\frac{y^{0,46} -1}{0,46}\), \(x_1\) = Orçamento em dólares, \(x_2\) = Quantidade de avaliações dos usuários.

Agora os regressores são significantes, e o modelo apresenta um \(R^2 \approx R^2_{adj} \approx 0,57\) indicando que cerca de 57% da variação da variável resposta é explicada pelos regressores. A interpretabilidade do modelo torna-se mais complicada devido a transformação realizada, mas temos que em média o incremento de uma unidade de \(x_1\) (mantendo \(x_2\) fixo) ocasiona um aumento de cerca de \(4,41.10^{-5}\) na resposta \(\frac{Bilheteria^{0,46} -1}{0,46}\). Enquanto que em média o incremento de uma unidade de \(x_2\) (mantendo \(x_1\) fixo) ocasiona um aumento de cerca de \(8,52.10^{-3}\) na resposta \(\frac{Bilheteria^{0,46} -1}{0,46}\).

Diagnóstico

Para todos os testes a seguir utilizaremos um nível de significância de \(\alpha = 0,05\). Primeiro realizamos o teste RESET de Ramsey que forneceu um p-valor superior a \(0,05\), portanto não rejeitamos a hipótese nula de linearidade do modelo, indicando que o modelo está corretamente especificado.

## 
##  RESET test
## 
## data:  ajuste_t2
## RESET = 2.9362, df1 = 2, df2 = 115, p-value = 0.05706

Através do Q-Q plot normal dos resíduos observamos que 7 pontos encontram-se fora dos envelopes, como \(\frac{7}{120} \approx 0,058\) então cerca de 5% das observações encontram-se fora dos envelopes, o que é razoável. Para garantir a normalidade dos resíduos realizamos o teste de Bera-Jarque, que forneceu um p-valor grande levando a não rejeição da hipótese nula de normalidade.

## Gaussian model (lm object)

## 
##  Jarque Bera Test
## 
## data:  res_t2
## X-squared = 3.8145, df = 2, p-value = 0.1485

A partir do gráfico de resíduos studentizados por índice a seguir verificamos um padrão aleatório, indicando a independência dos resíduos que tinhamos por suposição anteriormente. Já no gráfico de resíduos studentizados por valores ajustados também notamos certa aleatoriedade, indicando a homoscedasticidade.

Para testar a heteroscedasticidade utilizamos o teste de Goldfeld-Quandt, visto que sob normalidade ele costuma ser mais poderoso que os testes de Breusch-Pagan e Koenker. Omitindo \(\frac{1}{3}\) das observações centrais obtemos um p-valor grande indicando a não rejeição da hipótese nula de homoscedasticidade.

## 
##  Goldfeld-Quandt test
## 
## data:  ajuste_t2
## GQ = 1.3053, df1 = 37, df2 = 37, p-value = 0.2108
## alternative hypothesis: variance increases from segment 1 to 2

Apêndice: Código

library(tidyverse)
library(readr)
library(MASS)
library(psych)
library(tseries)
library(car)
library(lmtest)
library(hnp)
library(olsrr)

movie_metadata <- read_csv("movie_metadata.csv")

int <- movie_metadata

movie_metadata <- movie_metadata %>% filter(!is.na(budget))

movie_metadata <- movie_metadata %>% filter(budget > 0)

movie_metadata <- movie_metadata %>% filter(gross > 0)

movie_metadata <- movie_metadata %>% filter(title_year == 2014)

movie_metadata <- movie_metadata %>% filter(language == "English")

movie_metadata <- movie_metadata %>% filter(country == "USA")

dados <- na.omit(movie_metadata)

dados <- dados[,c(3,4,5,6,8,9,13,14,19,22,23,25,26,28)]

dados$content_rating <- as.factor(dados$content_rating)

ggplot() + geom_histogram(aes(x=dados$gross), color = "white", bins = 30) + ylab("Frequência") + xlab("Bilheteria em dólares")

## MODELO 1
ajuste <- lm(gross ~., data = dados)
summary(ajuste)
## MODELO 2
ajuste_AIC <- stepAIC(ajuste)


## MODELO 2
summary(ajuste_AIC)

dados <- dados[,c(3, 6, 7, 11, 14)]

res_AIC <- resid(ajuste_AIC)
fit_AIC <- fitted.values(ajuste_AIC)
stud.res_AIC <- studres(ajuste_AIC)

pairs.panels(x = dados, ellipses = FALSE, smooth = FALSE)

jarque.bera.test(res_AIC)
bptest(ajuste_AIC)

transfor <- boxcox(ajuste_AIC)
lambda <- transfor$x[which.max(transfor$y)]
lambda

bilhet <- dados$gross
bilhet_t <- bcPower(bilhet, lambda)

dados_t <- dados
dados_t$gross_t <- bilhet_t

ajuste_t <- lm(gross_t ~ director_facebook_likes + num_voted_users + budget + movie_facebook_likes, data = dados_t)
summary(ajuste_t)
res_t <- resid(ajuste_t)
fit_t <- fitted.values(ajuste_t)
stud.res_t <- studres(ajuste_t)
jarque.bera.test(res_t)

bptest(ajuste_t)

ajuste_t2 <- lm(gross_t ~ budget + num_voted_users, data = dados_t)
summary(ajuste_t2)
res_t2 <- resid(ajuste_t2)
fit_t2 <- fitted.values(ajuste_t2)
stud.res_t2 <- studres(ajuste_t2)

resettest(ajuste_t2)

hnp(ajuste_t2, halfnormal = F, resid.type = 'standard')
jarque.bera.test(res_t2)

dados_t %>% ggplot() + 
  geom_point(aes(x = 1:nrow(dados_t), y = stud.res_t2)) + 
  geom_abline(intercept = 2, slope = 0, color = 'red') +
  geom_abline(intercept = -2, slope = 0, color = 'red') + 
  xlab("Índice") + ylab("Resíduos Studentizados")

dados_t %>% ggplot() + 
  geom_point(aes(x = fit_t2, y = stud.res_t2)) +
  geom_abline(intercept = 2, slope = 0, color = 'red') +
  geom_abline(intercept = -2, slope = 0, color = 'red') +xlab("Valores ajustados") +ylab("Resíduos Studentizados")

gqtest(ajuste_t2, fraction = 1/3, order.by = dados_t$budget)