Load packages

library(MASS)
library(tidyverse)
library(statsr)
library(BAS)
library(broom)

Load data

load("movies.Rdata")

Part 1: Data

As it was mentioned first that the data set is comprised of 651 randomly sampled movies produced and released before 2016, we can conclude that the conclusion derived are generalizable and can be applied to the population.

It is not causal as there is no random assignment. So we can’t find reason on why people like a movie or not.

Part 2: Data manipulation

creating additional variables

movies <- movies %>% mutate(feature_film = ifelse(title_type == "Feature Film", "yes", "no"))
movies <- movies %>% mutate(drama = ifelse(genre == "Drama", "yes", "no"))
movies <- movies %>% mutate(mpaa_rating_R = ifelse(mpaa_rating == "R", "yes", "no"))
movies <- movies %>% mutate(oscar_season = ifelse(thtr_rel_month == 10 | thtr_rel_month == 11 | thtr_rel_month == 12, "yes", "no"))
movies <- movies %>% mutate(summer_season = ifelse(thtr_rel_month == 5 | thtr_rel_month == 6 | thtr_rel_month == 7 | thtr_rel_month == 8, "yes", "no"))

str(movies$feature_film)

##  chr [1:651] "yes" "yes" "yes" "yes" "yes" "no" "yes" "yes" "no" "yes" ...

str(movies$drama)

##  chr [1:651] "yes" "yes" "no" "yes" "no" "no" "yes" "yes" "no" "yes" ...

str(movies$mpaa_rating_R)

##  chr [1:651] "yes" "no" "yes" "no" "yes" "no" "no" "yes" "no" "no" ...

str(movies$oscar_season)

##  chr [1:651] "no" "no" "no" "yes" "no" "no" "no" "yes" "no" "no" "no" ...

str(movies$summer_season)

##  chr [1:651] "no" "no" "yes" "no" "no" "no" "no" "no" "no" "no" "yes" ...

Part 3: Exploratory data analysis

Part 3.1 : Between seasons and audience score

ggplot(movies, aes(x = factor(summer_season), y = audience_score)) + geom_boxplot()

ggplot(movies, aes(x = factor(oscar_season), y = audience_score)) + geom_boxplot()

Seems there is no difference based on the season in which the movie is released, as the plot looks similar in both the cases. So these variables are not at all related to each other.

Part 3.2 : Between feature_film and audience score

ggplot(movies, aes(x = factor(feature_film), y = audience_score)) + geom_boxplot()

Compared to the ‘feature films’ other types like ‘Documentaries and TV movies’ get much higher rating

Part 3.3 : Audience score vs drama

movies %>% group_by(drama) %>% summarise(min=min(audience_score), max=max(audience_score), iqr=IQR(audience_score), sd=sd(audience_score))

## # A tibble: 2 x 5
##   drama   min   max   iqr    sd
##   <chr> <dbl> <dbl> <dbl> <dbl>
## 1 no       11    97    38  21.3
## 2 yes      13    95    28  18.5

ggplot(movies, aes(x = factor(drama), y = audience_score)) + geom_boxplot()

In case of Drama, people preferred it a bit more than non-drama movies (based on 25-75% quantiles)

Part 3.4 : Audience score vs mpaa_rating_R

movies %>% group_by(mpaa_rating_R) %>% summarise(min=min(audience_score), max=max(audience_score), iqr=IQR(audience_score), sd=sd(audience_score))

## # A tibble: 2 x 5
##   mpaa_rating_R   min   max   iqr    sd
##   <chr>         <dbl> <dbl> <dbl> <dbl>
## 1 no               11    96  31.8  20.3
## 2 yes              14    97  35    20.2

ggplot(movies, aes(x = factor(mpaa_rating_R), y = audience_score)) + geom_boxplot()

Seems not much of a difference in the audience score

Part 4: Modeling

Part 4.1 : Model Selection

For this, we will use backward elimination with BIC . Let’s include all the params and get the BIC.

movies_no_na <- na.omit(movies)

mlm_movie <- lm(audience_score ~ feature_film + drama + runtime + mpaa_rating_R + thtr_rel_year + oscar_season + summer_season + imdb_rating + imdb_num_votes + critics_score + best_pic_nom + best_pic_win + best_actor_win + best_actress_win + best_dir_win + top200_box, data = movies_no_na)

BIC(mlm_movie)

## [1] 4716.774

BIC is 4716.774

BIC(lm(audience_score ~ feature_film + drama + runtime + mpaa_rating_R + thtr_rel_year + oscar_season + summer_season + imdb_rating + imdb_num_votes + critics_score + best_pic_nom + best_pic_win + best_actor_win + best_actress_win + best_dir_win, data = movies_no_na))

## [1] 4710.428

BIC is 4710.428. So remove the variable top200_box

BIC(lm(audience_score ~ feature_film + drama + runtime + mpaa_rating_R + thtr_rel_year + summer_season + imdb_rating + imdb_num_votes + critics_score + best_pic_nom + best_pic_win + best_actor_win + best_actress_win + best_dir_win, data = movies_no_na))

## [1] 4704.099

BIC is 4704.099, so remove the variable oscar_season

BIC(lm(audience_score ~ feature_film + drama + runtime + mpaa_rating_R + thtr_rel_year + summer_season + imdb_rating + imdb_num_votes + critics_score + best_pic_nom + best_actor_win + best_actress_win + best_dir_win, data = movies_no_na))

## [1] 4698.11

BIC is 4698.11, so remove the variable best_pic_win

BIC(lm(audience_score ~ feature_film + drama + runtime + mpaa_rating_R + thtr_rel_year + summer_season + imdb_rating + imdb_num_votes + critics_score + best_pic_nom + best_actor_win + best_actress_win, data = movies_no_na))

## [1] 4692.601

BIC is 4692.601, so remove the variable best_dir_win

BIC(lm(audience_score ~ feature_film + drama + runtime + mpaa_rating_R + thtr_rel_year + summer_season + imdb_rating + imdb_num_votes + critics_score + best_pic_nom + best_actress_win, data = movies_no_na))

## [1] 4687.432

BIC is 4687.432, so remove the variable best_actor_win

BIC(lm(audience_score ~ feature_film + drama + runtime + mpaa_rating_R + thtr_rel_year + imdb_rating + imdb_num_votes + critics_score + best_pic_nom + best_actress_win, data = movies_no_na))

## [1] 4682.776

BIC is 4682.776, so remove the variable summer_season

BIC(lm(audience_score ~ drama + runtime + mpaa_rating_R + thtr_rel_year + imdb_rating + imdb_num_votes + critics_score + best_pic_nom + best_actress_win, data = movies_no_na))

## [1] 4678.23

BIC is 4678.23, so remove the variable feature_film

BIC(lm(audience_score ~ runtime + mpaa_rating_R + thtr_rel_year + imdb_rating + imdb_num_votes + critics_score + best_pic_nom + best_actress_win, data = movies_no_na))

## [1] 4673.041

BIC is 4673.041, so remove the variable drama

BIC(lm(audience_score ~ runtime + mpaa_rating_R + thtr_rel_year + imdb_rating + critics_score + best_pic_nom + best_actress_win, data = movies_no_na))

## [1] 4668.281

BIC is 4668.281, so remove the variable imdb_num_votes

BIC(lm(audience_score ~ runtime + mpaa_rating_R + imdb_rating + critics_score + best_pic_nom + best_actress_win, data = movies_no_na))

## [1] 4663.531

BIC is 4663.531, so remove the variable thtr_rel_year

BIC(lm(audience_score ~ runtime + mpaa_rating_R + imdb_rating + critics_score + best_pic_nom, data = movies_no_na))

## [1] 4660.204

BIC is 4660.204, so remove the variable best_actress_win

BIC(lm(audience_score ~ runtime + mpaa_rating_R + imdb_rating + critics_score, data = movies_no_na))

## [1] 4656.76

BIC is 4656.76, so remove the variable best_pic_nom

BIC(lm(audience_score ~ runtime + imdb_rating + critics_score, data = movies_no_na))

## [1] 4653.581

BIC is 4653.581, so remove the variable mpaa_rating_R

We can’t remove anything else. Our final variables are runtime, critics_score and imdb_rating

Part 4.2 : Model diagnostics

The Bayesian model specification assumes that the errors are normally distributed with a constant variance and that the mean expected weekly wages is linear in IQ. We can check these assumption by examining the distribution of the residuals for the model.

In order to do so we will use predicted values, residuals, and standardized residuals of the model we fit earlier. The augment function in the broom package is going to come in handy here as it takes in a model object (the output of an lm) and returns a data frame with columns correspinding to variables in the model as well as predicted values (.fitted), residuals (.resid), and standardized residuals (.std.resid), along with a few others.

mlm_movie_model <- lm(audience_score ~ runtime + imdb_rating + critics_score, data = movies_no_na)
summary(mlm_movie_model)

## 
## Call:
## lm(formula = audience_score ~ runtime + imdb_rating + critics_score, 
##     data = movies_no_na)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.996  -6.706   0.738   5.435  52.550 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -32.90825    3.35451  -9.810  < 2e-16 ***
## runtime        -0.05791    0.02238  -2.588 0.009891 ** 
## imdb_rating    14.95043    0.60249  24.814  < 2e-16 ***
## critics_score   0.07531    0.02227   3.381 0.000769 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.15 on 615 degrees of freedom
## Multiple R-squared:  0.7483, Adjusted R-squared:  0.747 
## F-statistic: 609.4 on 3 and 615 DF,  p-value: < 2.2e-16

mlm_movie_model_aug <- augment(mlm_movie_model)

Part 4.2.1 : Linearity and Constant Variance:

We can verify this condition with a plot of the residuals vs. fitted (predicted) values.

ggplot(data = mlm_movie_model_aug, aes(x = .fitted, y = .resid)) +
  geom_point(alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Fitted values", y = "Residuals")

Part 4.2.2 : Normality:

To check this condition, we can look at a histogram of residuals

ggplot(data = mlm_movie_model_aug, aes(x = .resid)) +
  geom_histogram() +
  xlab("Residuals")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

or a normal probability plot of the residuals

ggplot(mlm_movie_model_aug) +
  geom_qq(aes(sample = .std.resid)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  labs(x = "Theoretical quantiles", y = "Standardized residuals")

Part 4.2.3 : Outliers

outliers <- Bayes.outlier(mlm_movie_model, k = 3)
outliers_df <- data.frame(probability = outliers$prob.outlier,
                          case = 1:length(outliers$prob.outlier))
ggplot(outliers_df, aes(ymax = probability, x = case)) +
  geom_linerange(ymin = 0) +
  labs(y = "Probability")

Outliers with prob > 0.5 is

outliers_df %>%
  filter(probability > 0.50)

## [1] probability case       
## <0 rows> (or 0-length row.names)

All conditions are checked. We are good to go.

Part 4.3 : Interpretation of the model co-efficients

Build model with the predicted params

bma_movie_model <- bas.lm(audience_score ~ runtime + imdb_rating + critics_score, data = movies_no_na,
                   prior = "BIC", 
                   modelprior = uniform())

bma_movie_model

## 
## Call:
## bas.lm(formula = audience_score ~ runtime + imdb_rating + critics_score, 
##     data = movies_no_na, prior = "BIC", modelprior = uniform())
## 
## 
##  Marginal Posterior Inclusion Probabilities: 
##     Intercept        runtime    imdb_rating  critics_score  
##        1.0000         0.5397         1.0000         0.9323

coef_movie <- coefficients(bma_movie_model)
coef_movie

## 
##  Marginal Posterior Summaries of Coefficients: 
## 
##  Using  BMA 
## 
##  Based on the top  8 models 
##                post mean  post SD   post p(B != 0)
## Intercept      62.21002    0.40900   1.00000      
## runtime        -0.03138    0.03333   0.53966      
## imdb_rating    14.89079    0.72722   1.00000      
## critics_score   0.07129    0.02889   0.93233

From the co-efficients, for the prediction of the audience_score we need intercept 100%, imdb_rating 100%, critics_score 93.233% and with a very less 53.966% runtime.

Same can be visually interpreted as

plot(coef_movie, subset = c(1,2,3,4), ask = FALSE)

Part 5: Prediction

Let’s predict the audience_score for the 2016 movie ROGUE ONE: A STAR WARS STORY using BMA

First check if it already there in the sample

movies %>% filter(grepl("rogue", title) | grepl("star", title)) %>% select(title)

## # A tibble: 1 x 1
##   title                  
##   <chr>                  
## 1 Bastard out of Carolina

No, it is safe to proceed.

movie_pred <- data.frame(runtime = c(133), imdb_rating=c(7.8), critics_score=c(84))
bma_predict <- predict(bma_movie_model, newdata=movie_pred, se.fit = TRUE, estimator = "BMA")
ci_bma_predict_movie <- confint(bma_predict, estimator = "BMA")
ci_bma_predict_movie

##          2.5%    97.5%     pred
## [1,] 63.15778 102.8948 82.84599
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"

Original Audience Score is 86%

Predicted Audience Score is 83%

It is a acceptable score. Hence our model works

References : IMDB Link and Rotten tomatoes link

Rotten Tomotoes IMDB

Part 6: Conclusion

All else held constants, having the imdb_rating (with 100%) in our model will increase the audience_score by 14.95043 compared to increase of only 0.07531 by critics_score (with 93.233%). So it is clear that imdb rating (given by users) are weighed more than critics score
Some good params like mpaa_rating_R, best_pic_nom, best_pic_win has no voice over the audience_score

So to conclude, audience_score even though it is directly related to people’s opinion as critic???s score/imdb_rating, but factors like runtime plays a major role but not best_pic_win/mpaa_rating_R.

But there are cases where the model fails if the critics like the movie too much and film type is of non-drama, the audience score???s could be low due to other factors like bad movie plot. As adj-r-squared value is 74.7%, there are still hefty of chances where the prediction will still go wrong.

Table of Contents