Modeling and prediction for movies
copied version
September 8, 2018
Modeling and prediction for movies Setup Load packages:
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(statsr)Load data
load("movies.Rdata")Part 1: Data The movies data were collected from 651 randomly sampled movies before 2016. This is an oberservational data and only association relationship can be concluded, we can not say any causality from this dataset.
Part 2: Research question In this analysis, we want to find out what attributes make a movie popular? Also, at the end of the report, we will build a prediction model and pick a movie from 2016 to validate our model predition.
Part 3: Exploratory data analysis Let’s review the variables we have here:
names(movies)## [1] "title" "title_type" "genre"
## [4] "runtime" "mpaa_rating" "studio"
## [7] "thtr_rel_year" "thtr_rel_month" "thtr_rel_day"
## [10] "dvd_rel_year" "dvd_rel_month" "dvd_rel_day"
## [13] "imdb_rating" "imdb_num_votes" "critics_rating"
## [16] "critics_score" "audience_rating" "audience_score"
## [19] "best_pic_nom" "best_pic_win" "best_actor_win"
## [22] "best_actress_win" "best_dir_win" "top200_box"
## [25] "director" "actor1" "actor2"
## [28] "actor3" "actor4" "actor5"
## [31] "imdb_url" "rt_url"Our objective is to find which variables can provide us better prediction of movie popularity. As we can see we have several scores from either critics or audience, how do we decide which variable to be our response variable? Let’s check the correlation coefficient for critics_score and audience_score.
movies %>% summarise(cor(critics_score,audience_score))## # A tibble: 1 x 1
## `cor(critics_score, audience_score)`
## <dbl>
## 1 0.704Since these two variables are highly correlated, it is ok to use either one and I choose audience_score as my response variable in this analysis.
I want to make the data cleaner and remove some variables that I won’t need in this analysis. Check the structure of the data.
str(movies)## Classes 'tbl_df', 'tbl' and 'data.frame': 651 obs. of 32 variables:
## $ title : chr "Filly Brown" "The Dish" "Waiting for Guffman" "The Age of Innocence" ...
## $ title_type : Factor w/ 3 levels "Documentary",..: 2 2 2 2 2 1 2 2 1 2 ...
## $ genre : Factor w/ 11 levels "Action & Adventure",..: 6 6 4 6 7 5 6 6 5 6 ...
## $ runtime : num 80 101 84 139 90 78 142 93 88 119 ...
## $ mpaa_rating : Factor w/ 6 levels "G","NC-17","PG",..: 5 4 5 3 5 6 4 5 6 6 ...
## $ studio : Factor w/ 211 levels "20th Century Fox",..: 91 202 167 34 13 163 147 118 88 84 ...
## $ thtr_rel_year : num 2013 2001 1996 1993 2004 ...
## $ thtr_rel_month : num 4 3 8 10 9 1 1 11 9 3 ...
## $ thtr_rel_day : num 19 14 21 1 10 15 1 8 7 2 ...
## $ dvd_rel_year : num 2013 2001 2001 2001 2005 ...
## $ dvd_rel_month : num 7 8 8 11 4 4 2 3 1 8 ...
## $ dvd_rel_day : num 30 28 21 6 19 20 18 2 21 14 ...
## $ imdb_rating : num 5.5 7.3 7.6 7.2 5.1 7.8 7.2 5.5 7.5 6.6 ...
## $ imdb_num_votes : int 899 12285 22381 35096 2386 333 5016 2272 880 12496 ...
## $ critics_rating : Factor w/ 3 levels "Certified Fresh",..: 3 1 1 1 3 2 3 3 2 1 ...
## $ critics_score : num 45 96 91 80 33 91 57 17 90 83 ...
## $ audience_rating : Factor w/ 2 levels "Spilled","Upright": 2 2 2 2 1 2 2 1 2 2 ...
## $ audience_score : num 73 81 91 76 27 86 76 47 89 66 ...
## $ best_pic_nom : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ best_pic_win : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ best_actor_win : Factor w/ 2 levels "no","yes": 1 1 1 2 1 1 1 2 1 1 ...
## $ best_actress_win: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ best_dir_win : Factor w/ 2 levels "no","yes": 1 1 1 2 1 1 1 1 1 1 ...
## $ top200_box : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ director : chr "Michael D. Olmos" "Rob Sitch" "Christopher Guest" "Martin Scorsese" ...
## $ actor1 : chr "Gina Rodriguez" "Sam Neill" "Christopher Guest" "Daniel Day-Lewis" ...
## $ actor2 : chr "Jenni Rivera" "Kevin Harrington" "Catherine O'Hara" "Michelle Pfeiffer" ...
## $ actor3 : chr "Lou Diamond Phillips" "Patrick Warburton" "Parker Posey" "Winona Ryder" ...
## $ actor4 : chr "Emilio Rivera" "Tom Long" "Eugene Levy" "Richard E. Grant" ...
## $ actor5 : chr "Joseph Julian Soria" "Genevieve Mooy" "Bob Balaban" "Alec McCowen" ...
## $ imdb_url : chr "http://www.imdb.com/title/tt1869425/" "http://www.imdb.com/title/tt0205873/" "http://www.imdb.com/title/tt0118111/" "http://www.imdb.com/title/tt0106226/" ...
## $ rt_url : chr "//www.rottentomatoes.com/m/filly_brown_2012/" "//www.rottentomatoes.com/m/dish/" "//www.rottentomatoes.com/m/waiting_for_guffman/" "//www.rottentomatoes.com/m/age_of_innocence/" ...Imagine when you have a brand new movie, what you can know right away are the title, genre, runtime, actor/actress, and director, and these are your explanatory variables. On the other hand, you do not know what is the critics score, audience rating, etc. at that particular moment. So my new dataset will contain only the explanatory variables and one response variable.
movies_1 <- movies %>% select(1:5,thtr_rel_year,thtr_rel_month,audience_score, best_actor_win,best_actress_win,best_dir_win)
names(movies_1)## [1] "title" "title_type" "genre"
## [4] "runtime" "mpaa_rating" "thtr_rel_year"
## [7] "thtr_rel_month" "audience_score" "best_actor_win"
## [10] "best_actress_win" "best_dir_win"First, I want to validate the data to see if there is any bias. Since the movies were collected randomly from each year, I calculate how many movies were released each year in the dataset.
How many movies were release each year?
yr <- movies_1 %>% group_by(thtr_rel_year) %>% summarise(count = n())
ggplot(yr,aes(x=as.factor(thtr_rel_year),y=count)) + geom_bar(stat = "identity") +
scale_x_discrete(breaks = seq(1970,2014,5)) + labs(x = "theater release year",y = "number of movies")
Even though more movies after 1995 were selected in this data, the difference is not that significant. Since we do not know total number of movies each year, it is reasonable to believe that more movies were released after 2000 and it is likely the reason behind the difference.
One more thing I want to explore is: Is there any certain month that the average scores are higher? Does not find anything interesting.
movies_1 %>% group_by(thtr_rel_month) %>% summarise(count = n(), avg_score = mean(audience_score))## # A tibble: 12 x 3
## thtr_rel_month count avg_score
## <dbl> <int> <dbl>
## 1 1 69 62.9
## 2 2 34 63.9
## 3 3 51 59.7
## 4 4 45 61.1
## 5 5 44 62.4
## 6 6 72 61.7
## 7 7 48 62.1
## 8 8 44 61.1
## 9 9 53 61.7
## 10 10 70 62.0
## 11 11 51 61.4
## 12 12 70 67.1Let’s check whether certain types of movies title_type have higher scores on average. Documentary type has much higher scores than other two. Notice that majority of movies are belong to Feature Film.
movies_1 %>% group_by(title_type) %>% summarise(count = n(), avg_score = mean(audience_score))## # A tibble: 3 x 3
## title_type count avg_score
## <fct> <int> <dbl>
## 1 Documentary 55 83.3
## 2 Feature Film 591 60.5
## 3 TV Movie 5 56.8Let’s see which genre gets higher scores on average. Still, Documentary has the highest avergae score.
movies_1 %>% group_by(genre) %>% summarise(count = n(), avg_score = mean(audience_score)) %>% arrange(desc(avg_score))## # A tibble: 11 x 3
## genre count avg_score
## <fct> <int> <dbl>
## 1 Documentary 52 82.8
## 2 Musical & Performing Arts 12 80.2
## 3 Other 16 66.7
## 4 Drama 305 65.3
## 5 Art House & International 14 64
## 6 Animation 9 62.4
## 7 Mystery & Suspense 59 55.9
## 8 Action & Adventure 65 53.8
## 9 Comedy 87 52.5
## 10 Science Fiction & Fantasy 9 50.9
## 11 Horror 23 45.8How about whether or not one of the main actors/actresses or director ever won a oscar? Here we see an interesting result. The top 2 highest average scores are coming from either “best actor and director” or “best actress and director”, but not all three. Notice that majority of movies filmed without best actor,actress, and director, indicating that there are still a lot of movie makers not getting Oscar award yet.
movies_1 %>% group_by(best_actor_win,best_actress_win,best_dir_win) %>% summarise(count = n(), avg_score = mean(audience_score)) %>% arrange(desc(avg_score))## # A tibble: 8 x 5
## # Groups: best_actor_win, best_actress_win [4]
## best_actor_win best_actress_win best_dir_win count avg_score
## <fct> <fct> <fct> <int> <dbl>
## 1 yes no yes 8 76.9
## 2 no yes yes 5 72.2
## 3 yes yes yes 4 72
## 4 yes yes no 16 66.4
## 5 no no yes 26 66.3
## 6 no no no 480 62.0
## 7 no yes no 47 61.4
## 8 yes no no 65 60.3Part 4: Modeling As above mentioned, I removed some variables and only included 10 explanatory variables(known) and 1 response variable(unknown). Here I would like to use backwards elimination p-value methods since it’s the fastest method. In my full model, I exclude title and thtr_rel_year. It might be possible that people enjoy movies more when certain “words” appear in the title, but we don’t have skills to analyze at this point. For the release year, it’s not related when we predict certain movie is popular or not, so I didn’t include thtr_rell_year in my model.
# Full model
mlr <- lm(audience_score ~ title_type + genre + runtime + mpaa_rating + thtr_rel_month +
best_actor_win + best_actress_win + best_dir_win, data = movies_1)
summary(mlr)##
## Call:
## lm(formula = audience_score ~ title_type + genre + runtime +
## mpaa_rating + thtr_rel_month + best_actor_win + best_actress_win +
## best_dir_win, data = movies_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52.684 -12.995 1.278 13.251 39.492
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.05485 9.00873 6.444 2.32e-10 ***
## title_typeFeature Film -11.83293 6.61174 -1.790 0.07399 .
## title_typeTV Movie -21.40435 10.41533 -2.055 0.04028 *
## genreAnimation 4.67673 6.97068 0.671 0.50252
## genreArt House & International 9.53012 5.39195 1.767 0.07764 .
## genreComedy 0.98293 2.98616 0.329 0.74214
## genreDocumentary 16.62940 7.09879 2.343 0.01946 *
## genreDrama 10.82022 2.53592 4.267 2.29e-05 ***
## genreHorror -6.85914 4.44890 -1.542 0.12364
## genreMusical & Performing Arts 19.77580 6.05167 3.268 0.00114 **
## genreMystery & Suspense 1.07470 3.33231 0.323 0.74717
## genreOther 11.97678 5.04504 2.374 0.01790 *
## genreScience Fiction & Fantasy -4.02326 6.35636 -0.633 0.52700
## runtime 0.17768 0.04245 4.185 3.26e-05 ***
## mpaa_ratingNC-17 -10.41867 13.53800 -0.770 0.44183
## mpaa_ratingPG -9.76940 4.90043 -1.994 0.04663 *
## mpaa_ratingPG-13 -15.63729 5.00355 -3.125 0.00186 **
## mpaa_ratingR -9.99881 4.85316 -2.060 0.03979 *
## mpaa_ratingUnrated -6.48384 5.60508 -1.157 0.24780
## thtr_rel_month -0.07782 0.20458 -0.380 0.70379
## best_actor_winyes -1.00440 2.11506 -0.475 0.63504
## best_actress_winyes -0.11666 2.32657 -0.050 0.96002
## best_dir_winyes 6.15676 2.92599 2.104 0.03576 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.83 on 627 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2495, Adjusted R-squared: 0.2231
## F-statistic: 9.474 on 22 and 627 DF, p-value: < 2.2e-16Model selection Remove best_actress_win since the p-value is the highest one at 0.96002
mlr <- lm(audience_score ~ title_type + genre + runtime + mpaa_rating + thtr_rel_month +
best_actor_win + best_dir_win, data = movies_1)
summary(mlr)##
## Call:
## lm(formula = audience_score ~ title_type + genre + runtime +
## mpaa_rating + thtr_rel_month + best_actor_win + best_dir_win,
## data = movies_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52.66 -12.99 1.23 13.26 39.50
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.08741 8.97817 6.470 1.98e-10 ***
## title_typeFeature Film -11.83355 6.60648 -1.791 0.07374 .
## title_typeTV Movie -21.41802 10.40349 -2.059 0.03993 *
## genreAnimation 4.65993 6.95709 0.670 0.50323
## genreArt House & International 9.51736 5.38166 1.768 0.07747 .
## genreComedy 0.96844 2.96978 0.326 0.74446
## genreDocumentary 16.62061 7.09099 2.344 0.01939 *
## genreDrama 10.80326 2.51126 4.302 1.96e-05 ***
## genreHorror -6.86630 4.44308 -1.545 0.12275
## genreMusical & Performing Arts 19.77684 6.04683 3.271 0.00113 **
## genreMystery & Suspense 1.05620 3.30919 0.319 0.74970
## genreOther 11.96525 5.03579 2.376 0.01780 *
## genreScience Fiction & Fantasy -4.02515 6.35120 -0.634 0.52647
## runtime 0.17738 0.04199 4.224 2.76e-05 ***
## mpaa_ratingNC-17 -10.39950 13.52185 -0.769 0.44213
## mpaa_ratingPG -9.77016 4.89651 -1.995 0.04644 *
## mpaa_ratingPG-13 -15.63761 4.99957 -3.128 0.00184 **
## mpaa_ratingR -9.99467 4.84861 -2.061 0.03968 *
## mpaa_ratingUnrated -6.47792 5.59938 -1.157 0.24775
## thtr_rel_month -0.07793 0.20441 -0.381 0.70316
## best_actor_winyes -1.01135 2.10883 -0.480 0.63169
## best_dir_winyes 6.15251 2.92244 2.105 0.03566 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.82 on 628 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2495, Adjusted R-squared: 0.2244
## F-statistic: 9.941 on 21 and 628 DF, p-value: < 2.2e-16Remove thtr_rel_month since the p-value is the highest one at 0.70316
mlr <- lm(audience_score ~ title_type + genre + runtime + mpaa_rating +
best_actor_win + best_dir_win, data = movies_1)
summary(mlr)##
## Call:
## lm(formula = audience_score ~ title_type + genre + runtime +
## mpaa_rating + best_actor_win + best_dir_win, data = movies_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52.902 -13.126 1.301 13.301 39.384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 57.90767 8.95968 6.463 2.06e-10 ***
## title_typeFeature Film -11.81662 6.60184 -1.790 0.07395 .
## title_typeTV Movie -21.29456 10.39138 -2.049 0.04085 *
## genreAnimation 4.56065 6.94749 0.656 0.51178
## genreArt House & International 9.51314 5.37800 1.769 0.07739 .
## genreComedy 0.92043 2.96510 0.310 0.75634
## genreDocumentary 16.62488 7.08616 2.346 0.01928 *
## genreDrama 10.81328 2.50942 4.309 1.90e-05 ***
## genreHorror -6.87936 4.43992 -1.549 0.12178
## genreMusical & Performing Arts 19.75002 6.04231 3.269 0.00114 **
## genreMystery & Suspense 1.12324 3.30227 0.340 0.73386
## genreOther 12.04755 5.02774 2.396 0.01686 *
## genreScience Fiction & Fantasy -4.00936 6.34675 -0.632 0.52780
## runtime 0.17394 0.04099 4.244 2.53e-05 ***
## mpaa_ratingNC-17 -10.27603 13.50878 -0.761 0.44713
## mpaa_ratingPG -9.78875 4.89294 -2.001 0.04587 *
## mpaa_ratingPG-13 -15.58932 4.99457 -3.121 0.00188 **
## mpaa_ratingR -10.00672 4.84521 -2.065 0.03931 *
## mpaa_ratingUnrated -6.44191 5.59478 -1.151 0.25000
## best_actor_winyes -1.04551 2.10550 -0.497 0.61967
## best_dir_winyes 6.13159 2.91994 2.100 0.03613 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.81 on 629 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2493, Adjusted R-squared: 0.2254
## F-statistic: 10.44 on 20 and 629 DF, p-value: < 2.2e-16Remove best_actor_win since the p-value is the highest one at 0.61967
best_mlr <- lm(audience_score ~ title_type + genre + runtime + mpaa_rating +
best_dir_win, data = movies_1)
summary(best_mlr)##
## Call:
## lm(formula = audience_score ~ title_type + genre + runtime +
## mpaa_rating + best_dir_win, data = movies_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52.671 -13.017 1.343 13.316 39.399
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58.3175 8.9163 6.541 1.27e-10 ***
## title_typeFeature Film -11.8562 6.5974 -1.797 0.07280 .
## title_typeTV Movie -21.2172 10.3840 -2.043 0.04144 *
## genreAnimation 4.4799 6.9414 0.645 0.51891
## genreArt House & International 9.5847 5.3728 1.784 0.07492 .
## genreComedy 0.9012 2.9631 0.304 0.76111
## genreDocumentary 16.5874 7.0815 2.342 0.01947 *
## genreDrama 10.7557 2.5052 4.293 2.04e-05 ***
## genreHorror -6.8424 4.4366 -1.542 0.12352
## genreMusical & Performing Arts 19.7933 6.0381 3.278 0.00110 **
## genreMystery & Suspense 0.9625 3.2844 0.293 0.76958
## genreOther 11.9933 5.0236 2.387 0.01726 *
## genreScience Fiction & Fantasy -3.9078 6.3397 -0.616 0.53785
## runtime 0.1696 0.0400 4.239 2.58e-05 ***
## mpaa_ratingNC-17 -10.6730 13.4770 -0.792 0.42869
## mpaa_ratingPG -9.8586 4.8880 -2.017 0.04413 *
## mpaa_ratingPG-13 -15.6005 4.9915 -3.125 0.00186 **
## mpaa_ratingR -10.0048 4.8423 -2.066 0.03923 *
## mpaa_ratingUnrated -6.4048 5.5909 -1.146 0.25241
## best_dir_winyes 6.0833 2.9166 2.086 0.03740 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.8 on 630 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.249, Adjusted R-squared: 0.2264
## F-statistic: 10.99 on 19 and 630 DF, p-value: < 2.2e-16Try to remove title_type, genre, mpaa_rating one at a time since the slope of some levels are not significant.(I use 95% confident level). However, adjusted R square are all decreasing. Hence I got the best model already (adjusted R sqare: 0.2264).
Model diagnostics
# histogram of residuals...nearly normal residuals --- condition met
hist(best_mlr$residuals)
# qqnorm
qqnorm(best_mlr$residuals)
# residual plot. residuals vs. predicted...constant variability of residuals --- condition met
plot(best_mlr$residuals ~ mlr$fitted)
# independent residuals...there is not time series structure
plot(best_mlr$residuals)
Our best model has p-value: < 2.2e-16, which is smaller than 0.05. So we can say our model has at least one of the beta that is non-zero and the model as a whole is significant. We have adjusted R sqare: 0.2264, which tells us 22.6% of the variability of the data can be explained by our model.
The model tells us that if the movie was filmed by an Oscar winning director, its audience score will be higher 6 points on average when all else held constand. If the movie time increases by 1 minute, the audience score will be up 0.17 points when all else held constand.
For the title_type, genre, mpaa_rating variables, I’m not sure how to interpret these findings. Some levels of variables are not significant, so I’m not sure if we can interpret the same way as other significant variables.
Part 5: Prediction Pick a movie from 2016: The movie I pick is “La La Land”. title_type: Feature Film Rating: PG-13 Runtime: 128 minutes Genre: Comedy, Drama, Musical & Performing Arts.let’s use Musical & Performing Arts in our analysis Director: Damien Chazelle
lalaland <- data.frame(title_type = "Feature Film",genre = "Musical & Performing Arts",runtime = 128,mpaa_rating = "PG-13",best_dir_win = "yes")
predict(best_mlr,lalaland)## 1
## 78.44126predict(best_mlr,lalaland,interval = "prediction",level = 0.95)## fit lwr upr
## 1 78.44126 41.33339 115.5491Our model predicts La La Land’s audience score to be 72.36. We have 95% confidence that the score will be within 35.59 and 109.13. Actual audience score: 81
Reference: Wikipedia: https://en.wikipedia.org/wiki/La_La_Land_(film) Rotten tomatoes: https://www.rottentomatoes.com/m/la_la_land/ * *
Part 6: Conclusion Overall, our model provides a way to predict audience score of a particular movie, based on title_type,genre,runtime,mpaa_rating,best_dir_win. It tells us if a movie were filmed by an award winning director, the score tends to be higher. However, some features of this model can not be explained very well, for example, why the higher movie runtime gets higher score? Also, some levels of variables are not significant in the model. I am not sure if we should keep them or not.