Modeling and prediction for movies
Setup
Load packages
library(ggplot2)
library(dplyr)
library(statsr)
library(grid)
library(gridExtra)
library(knitr)Load data
Make sure your data and R Markdown files are in the same directory. When loaded your data file will be called movies. Delete this note when before you submit your work.
load("movies.Rdata")Part 1: Data
The data set is comprised of 651 randomly sampled movies produced and released before 2016.
Some of these variables are only there for informational purposes and do not make any sense to include in a statistical analysis. It is up to you to decide which variables are meaningful and which should be omitted. For example information in the the actor1 through actor5 variables was used to determine whether the movie casts an actor or actress who won a best actor or actress Oscar.
You might also choose to omit certain observations or restructure some of the variables to make them suitable for answering your research questions.
When you are fitting a model you should also be careful about collinearity, as some of these variables may be dependent on each other. A possible non-independent bias may arise with regard to movie sequels (for example, all of the Star Wars episodes) whereby the popularity of a sequel movie may be influenced by that of the prior releases. A 2011 study by metacritic showed that critics and audiences tended to rate sequels and other derivative works higher than the original work. On the other hand, a 2016 chart produced by reddit of ratings for original movies versus their sequls (considering 1990-2016 films only) appears to indicate just the opposite. In light of
these conflicting results, not knowing how sequels were treated in the data sampling method, and not having a definitive means of identifying sequels and originals in the raw data set so as to try to identify a bias, it is assumed that no bias exists for the purpose of this analysis.
Overall, it is assumed that the individual movie observations are independent of each other. The size of the general population of movies produced is difficult to quantify exactly, however the IMDB web site lists over 105,000 feature film titles produced in the U.S. alone. A sample of 651 movies is certainly well below the 10% of the population threshold needed for independence.
As the raw data is a sample taken from existing data, this is an observational study and no causal relationships can be inferred or assumed from the conclusions drawn.
Part 2: Research question
The research question addressed in this analysis is this: is it possible to predict the popularity of a feature film prior to its release based upon certain characteristic information about it? More specifically, do variables such as movie genre, MPAA rating, run length, etc. work as reasonable predictors of a movie’s popularity?
A good predictive model would be valuable to a theater owner when choosing which movies to show. Popular movies are more likely to draw larger crowds of patrons, leading to higher ticket and concession revenues. Predicting which movies are likely to be popular ahead of their release would be valuable information when deciding which movies to show.
Part 3: Exploratory data analysis
Codebook: title: Title of movie title_type: Type of movie (Documentary, Feature Film, TV Movie) genre: Genre of movie (Action & Adventure, Comedy, Documentary, Drama, Horror, Mystery & Suspense, Other) runtime: Runtime of movie (in minutes) mpaa_rating: MPAA rating of the movie (G, PG, PG-13, R, Unrated) studio: Studio that produced the movie thtr_rel_year: Year the movie is released in theaters thtr_rel_month: Month the movie is released in theaters thtr_rel_day: Day of the month the movie is released in theaters dvd_rel_year: Year the movie is released on DVD dvd_rel_month: Month the movie is released on DVD dvd_rel_day: Day of the month the movie is released on DVD imdb_rating: Rating on IMDB imdb_num_votes: Number of votes on IMDB critics_rating: Categorical variable for critics rating on Rotten Tomatoes (Certified Fresh, Fresh, Rotten) critics_score: Critics score on Rotten Tomatoes audience_rating: Categorical variable for audience rating on Rotten Tomatoes (Spilled, Upright) audience_score: Audience score on Rotten Tomatoes best_pic_nom: Whether or not the movie was nominated for a best picture Oscar (no, yes) best_pic_win: Whether or not the movie won a best picture Oscar (no, yes) best_actor_win: Whether or not one of the main actors in the movie ever won an Oscar (no, yes) - note that this is not necessarily whether the actor won an Oscar for their role in the given movie best_actress win: Whether or not one of the main actresses in the movie ever won an Oscar (no, yes) - not that this is not necessarily whether the actresses won an Oscar for their role in the given movie best_dir_win: Whether or not the director of the movie ever won an Oscar (no, yes) - not that this is not necessarily whether the director won an Oscar for the given movie top200_box: Whether or not the movie is in the Top 200 Box Office list on BoxOfficeMojo (no, yes) director: Director of the movie actor1: First main actor/actress in the abridged cast of the movie actor2: Second main actor/actress in the abridged cast of the movie actor3: Third main actor/actress in the abridged cast of the movie actor4: Fourth main actor/actress in the abridged cast of the movie actor5: Fifth main actor/actress in the abridged cast of the movie imdb_url: Link to IMDB page for the movie rt_url: Link to Rotten Tomatoes page for the movie
There are 651 movies represented in the raw data. The following charts show a breakdown of the type of movies included in the sample.
p1 <- ggplot(data=movies, aes(x=genre)) +
geom_bar(fill="blue") +
xlab("Movie Genre") +
theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0))
p2 <- ggplot(data=movies, aes(x=title_type)) +
geom_bar(fill="blue") +
xlab("Movie Type") +
theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0))
p3 <- ggplot(data=movies, aes(x=mpaa_rating)) +
geom_bar(fill="blue") +
xlab("Movie MPAA Rating") +
theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0))
p4 <- ggplot(data=movies, aes(x=runtime)) +
geom_histogram(binwidth=10, fill="blue") +
xlab("Movie Runtime")
grid.arrange(p2, p3, p1, p4, nrow=2,
top="Movie Characteristics")## Warning: Removed 1 rows containing non-finite values (stat_bin).
There are 60 movies in the raw data that are type “Documentary” or “TV Movie”. These will be removed from the analysis as they would not likely be shown in a typical movie theater. Also, there are 52 movies with MPAA ratings of NC-17 or are unrated. These, as well, would not likely be shown in a typical movie theater and will be excluded from the analysis. Clean the data to remove the movies to be excluded from the analysis.
newmovies <- movies %>% filter(title_type=="Feature Film") %>%
filter(!(mpaa_rating %in% c("NC-17", "Unrated")))After removing these data, 573 movies remain in the processed data set for analysis. * * *
Part 4: Modeling
The target response variable for the prediction model is a movie rating score, but with 3 to choose from, which one should be used? Two of the ratings come from the Rotten Tomatoes web site: one is an average of reviews by movie critics and the other is an average of reviews from the audience. The third rating is an average of reviews on the IMDB web site. Following plot shows correlation between the different rating scores.
Note:Helper function for adding correlation coeficient values to a pairwise plot
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
# Using the helper function,the linear correlation R value is included in chart
pairs(~ imdb_rating + critics_score + audience_score,
data=newmovies, lower.panel=panel.smooth, upper.panel=panel.cor)
So,we use one rating as the response variable.
Histograms of the data for the three ratings:
p1 <- ggplot(data=newmovies, aes(x=imdb_rating)) +
geom_histogram(binwidth=0.5, fill="blue") +
xlab("IMDB Scores")
p2 <- ggplot(data=newmovies, aes(x=critics_score)) +
geom_histogram(binwidth=5, fill="blue") +
xlab("Critics Scores")
p3 <- ggplot(data=newmovies, aes(x=audience_score)) +
geom_histogram(binwidth=5, fill="blue") +
xlab("Audience Scores")
grid.arrange(p1, p2, p3, nrow=1,
top="Distribution of Rating Scores")
The IMDB scores show a nice, mostly normal distribution centered around a mean of 6.37 with somewhat of a left-skewed. Given its distribution and the fact that it has the highest pairwise correlation with the other scores, the IMDB rating was the chosen response variable.
Since the goal is to predict the popularity of a movie prior to its release, the prediction model uses only variables from the data set that could be known ahead of time. Thus, variables such as DVD release date, number of IMDB votes, best picture nomination/win, etc. were not chosen to be in the model. Variables with large domains, such as studio name, actor/director names, URLs, etc. were excluded as well.
The first pass model was a multi-variable linear regression model using all relevant variables as predictors to derive a predicted IMDB rating score. Specifically, the initial set of predictors was: genre runtime mpaa_rating thtr_rel_month best_actor_win best_actress win best_dir_win Theater release month was included assuming that movies released at certain times of the year may be more popular than others. Release year was discarded as being irrelevant and release day was thought to be too granular to be a worthwhile predictor.
Model fitting results:
modelfit <- lm(imdb_rating ~ genre + runtime + mpaa_rating + thtr_rel_month +
best_actor_win + best_actress_win + best_dir_win, data=newmovies)
summary(modelfit)##
## Call:
## lm(formula = imdb_rating ~ genre + runtime + mpaa_rating + thtr_rel_month +
## best_actor_win + best_actress_win + best_dir_win, data = newmovies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8670 -0.5307 0.0430 0.6227 1.9998
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.033861 0.376917 13.355 < 2e-16 ***
## genreAnimation -0.306365 0.373898 -0.819 0.412921
## genreArt House & International 0.386649 0.321918 1.201 0.230233
## genreComedy -0.071803 0.158142 -0.454 0.649977
## genreDocumentary 0.790707 0.675727 1.170 0.242441
## genreDrama 0.593064 0.134232 4.418 1.20e-05 ***
## genreHorror -0.113862 0.241935 -0.471 0.638089
## genreMusical & Performing Arts 0.922245 0.353553 2.609 0.009339 **
## genreMystery & Suspense 0.384125 0.175588 2.188 0.029113 *
## genreOther 0.736434 0.279093 2.639 0.008558 **
## genreScience Fiction & Fantasy -0.266850 0.333834 -0.799 0.424431
## runtime 0.015921 0.002688 5.923 5.56e-09 ***
## mpaa_ratingPG -0.790634 0.284121 -2.783 0.005574 **
## mpaa_ratingPG-13 -1.057501 0.289001 -3.659 0.000277 ***
## mpaa_ratingR -0.715151 0.282431 -2.532 0.011613 *
## thtr_rel_month 0.006145 0.011526 0.533 0.594129
## best_actor_winyes -0.044683 0.114659 -0.390 0.696908
## best_actress_winyes 0.065389 0.124816 0.524 0.600569
## best_dir_winyes 0.358104 0.155371 2.305 0.021545 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9365 on 554 degrees of freedom
## Multiple R-squared: 0.248, Adjusted R-squared: 0.2236
## F-statistic: 10.15 on 18 and 554 DF, p-value: < 2.2e-16anova(modelfit)## Analysis of Variance Table
##
## Response: imdb_rating
## Df Sum Sq Mean Sq F value Pr(>F)
## genre 10 95.24 9.524 10.8585 < 2.2e-16 ***
## runtime 1 41.21 41.211 46.9851 1.913e-11 ***
## mpaa_rating 3 18.54 6.181 7.0468 0.0001178 ***
## thtr_rel_month 1 0.25 0.252 0.2869 0.5924429
## best_actor_win 1 0.10 0.096 0.1093 0.7410594
## best_actress_win 1 0.27 0.272 0.3101 0.5778658
## best_dir_win 1 4.66 4.659 5.3122 0.0215453 *
## Residuals 554 485.92 0.877
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1With an adjusted R-square value of 0.2236, there is plenty of room for improvement in the model. Next is to one-by-one remove insignificant variables based on p values, eliminating the variable with the highest p value each time, until all remaining variables were significant.
P value was used for the elimination criteria because the goal was to create a model with the highest predictive value using only variables with demonstrable sigificance. It seems like a possible case of overfitting to create a model with an artificially inflated R-squared value due to incuding statistically insignificant predictors. The model results are summarized below:
modelfit2 <- lm(imdb_rating ~ genre + runtime + mpaa_rating +
best_dir_win, data=newmovies)
summary(modelfit2)##
## Call:
## lm(formula = imdb_rating ~ genre + runtime + mpaa_rating + best_dir_win,
## data = newmovies)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8689 -0.5436 0.0404 0.6056 2.0432
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.035253 0.368299 13.672 < 2e-16 ***
## genreAnimation -0.290502 0.371977 -0.781 0.435153
## genreArt House & International 0.392689 0.320509 1.225 0.221017
## genreComedy -0.059992 0.156688 -0.383 0.701958
## genreDocumentary 0.779294 0.673510 1.157 0.247742
## genreDrama 0.599714 0.132672 4.520 7.55e-06 ***
## genreHorror -0.107019 0.241248 -0.444 0.657499
## genreMusical & Performing Arts 0.924963 0.352691 2.623 0.008965 **
## genreMystery & Suspense 0.382923 0.173083 2.212 0.027347 *
## genreOther 0.731349 0.277479 2.636 0.008631 **
## genreScience Fiction & Fantasy -0.263124 0.332956 -0.790 0.429709
## runtime 0.016264 0.002488 6.537 1.42e-10 ***
## mpaa_ratingPG -0.790039 0.283350 -2.788 0.005481 **
## mpaa_ratingPG-13 -1.060774 0.288355 -3.679 0.000257 ***
## mpaa_ratingR -0.715358 0.281770 -2.539 0.011394 *
## best_dir_winyes 0.359052 0.155015 2.316 0.020907 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9346 on 557 degrees of freedom
## Multiple R-squared: 0.2471, Adjusted R-squared: 0.2268
## F-statistic: 12.19 on 15 and 557 DF, p-value: < 2.2e-16anova(modelfit2)## Analysis of Variance Table
##
## Response: imdb_rating
## Df Sum Sq Mean Sq F value Pr(>F)
## genre 10 95.24 9.524 10.9040 < 2.2e-16 ***
## runtime 1 41.21 41.211 47.1820 1.735e-11 ***
## mpaa_rating 3 18.54 6.181 7.0763 0.0001129 ***
## best_dir_win 1 4.69 4.686 5.3650 0.0209069 *
## Residuals 557 486.51 0.873
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Note that the adjusted R-square value (0.2268) is virtually the same as the initial model, but all predictor variables in the final model are statistically significant. With a virtually zero overall p value, we would reject a null hpothesis that all of the model coefficients are actually zero in favor of an alternative hypothesis that at least one of the coefficients is not.
The model raises the predicted movie rating if the director has ever won an Oscar (all other predictors held constant) whereas the variables for the same being true of the lead actor or actress were removed from the model as being statistically insignificant.
Model Diagnostics
pMod <- fortify(modelfit2)
# Create residuals scatter plot.
p1 <- ggplot(pMod, aes(x=.fitted, y=.resid))+geom_point() +
geom_smooth(se=FALSE)+geom_hline(yintercept=0, col="red", linetype="dashed") +
xlab("Fitted Values")+ylab("Residuals") +
ggtitle("Residual vs Fitted Plot")
pMod$.qqnorm <- qqnorm(pMod$.stdresid, plot.it=FALSE)$x
y <- quantile(pMod$.stdresid, c(0.25, 0.75)) # Find the 1st and 3rd quartiles
x <- quantile(pMod$.qqnorm, c(0.25, 0.75)) # Find the 1st and 3rd quartiles
slope <- diff(y) / diff(x) # Compute the line slope
int <- y[1] - slope * x[1] # Compute the line intercept
# Residuals QQ plot.
p2 <- ggplot(pMod, aes(.qqnorm, .stdresid)) +
geom_point(na.rm = TRUE) +
geom_abline(intercept=int, slope=slope, color="red") +
xlab("Theoretical Quantiles")+ylab("Standardized Residuals") +
ggtitle("Normal Q-Q Plot")
# Residuals histogram plot
p3 <- ggplot(data=pMod, aes(x=.resid)) +
geom_histogram(binwidth=0.5, fill="blue") +
xlab("Residuals") +
ggtitle("Distribution of Residuals")
grid.arrange(p1, p3, p2, nrow=1, top="Model Diagnostic Plots")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Part 5: Prediction
dataDG <- data.frame(genre="Comedy", runtime=102, mpaa_rating="R", best_dir_win="no")
predDG <- predict(modelfit2, dataDG, interval="predict")
dataDead <- data.frame(genre="Action & Adventure", runtime=108, mpaa_rating="R", best_dir_win="no")
predDead <- predict(modelfit2, dataDead, interval="predict")
# Show prediction results.
df <- data.frame(t=c("American Pie", "Ant-Man and the Wasp"),
p=c(sprintf("%2.1f", predDG[1]),
sprintf("%2.1f", predDead[1])),
i=c(sprintf("%2.1f - %2.1f", predDG[2], predDG[3]),
sprintf("%2.1f - %2.1f", predDead[2], predDead[3])),
r=c("6.0", "8.1"))
kable(df, col.names=c("Movie Title", "Predicted Rating", "95% Prediction Interval", "Actual Rating"))| Movie Title | Predicted Rating | 95% Prediction Interval | Actual Rating |
|---|---|---|---|
| American Pie | 5.9 | 4.1 - 7.8 | 6.0 |
| Ant-Man and the Wasp | 6.1 | 4.2 - 7.9 | 8.1 |
Show prediction results.
Note that the 95% confidence prediction intervals are very wide. This is a reflection of the poor predictive capability of the model (further evidenced by its F-statistic and adjusted R-square values). * * *
Part 6: Conclusion
In answer to the research question stated in Part 2, yes, it is possible to predict the popularity of a movie based upon basic movie characteristic data. A multi-variable linear regression model was used in prediction.
But, there is much room for improvement. As shown in Part 5, the predictive power of the model is limited.