Week 11 Data Dive
Jenna Horrall
2024-04-13
Week 11 Data Dive
Creating the Model
To begin creating a model, I first decided I wanted to examine gross sales. I currently do not know if a linear model or a gneralized linear model would be good for this scenario, so I just stared with a regular linear model. For my explanatory variables, I thought that the average rating of the book, sales rank, price of the book, and author rating could all be factors that influence how much money the book makes.
lm_model <- lm(`gross sales` ~ Book_average_rating + `sales rank` + `sale price` + Author_Rating, data = book_sales)
summary(lm_model)##
## Call:
## lm(formula = `gross sales` ~ Book_average_rating + `sales rank` +
## `sale price` + Author_Rating, data = book_sales)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5292 -1571 -676 748 40025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1883.0616 2807.1134 -0.671 0.50248
## Book_average_rating 1059.5598 653.0375 1.623 0.10499
## `sales rank` -4.5508 0.3035 -14.993 < 2e-16 ***
## `sale price` 338.5657 28.5028 11.878 < 2e-16 ***
## Author_RatingFamous -138.1847 515.1082 -0.268 0.78855
## Author_RatingIntermediate 933.9305 301.2118 3.101 0.00198 **
## Author_RatingNovice 3217.6443 822.3647 3.913 9.71e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3298 on 1063 degrees of freedom
## Multiple R-squared: 0.3024, Adjusted R-squared: 0.2984
## F-statistic: 76.79 on 6 and 1063 DF, p-value: < 2.2e-16ggplot(book_sales, aes(x = Book_average_rating, y = `gross sales`)) +
geom_point(color = "black") +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(x = "Book Average Rating", y = "Gross Sales") +
ggtitle("Scatterplot with Linear Regression Line") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
Based on the summary of the linear model and the scatter plot of the average book rating and gross sales, a linear model does not seem like the best fit in this scenario. Based on the R-squared value of .3024, the linear model is only capable of modeling about 30% of the variance of gross sales. For these reasons, I will continue with a generalized linear model.
Since gross sales is a form of count, I chose to indicate a poisson distribution.
str(book_sales)## spc_tbl_ [1,070 × 15] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ index : num [1:1070] 0 1 2 3 4 5 6 7 8 9 ...
## $ Publishing Year : num [1:1070] 1975 1987 2015 2008 2011 ...
## $ Book Name : chr [1:1070] "Beowulf" "Batman: Year One" "Go Set a Watchman" "When You Are Engulfed in Flames" ...
## $ Author : chr [1:1070] "Unknown, Seamus Heaney" "Frank Miller, David Mazzucchelli, Richmond Lewis, Dennis O'Neil" "Harper Lee" "David Sedaris" ...
## $ language_code : chr [1:1070] "en-US" "eng" "eng" "en-US" ...
## $ Author_Rating : chr [1:1070] "Novice" "Intermediate" "Novice" "Intermediate" ...
## $ Book_average_rating: num [1:1070] 3.42 4.23 3.31 4.04 4.04 4.08 4.03 3.9 4.34 4.2 ...
## $ Book_ratings_count : num [1:1070] 155903 145267 138669 150898 198283 ...
## $ genre : chr [1:1070] "fiction" "fiction" "fiction" "fiction" ...
## $ gross sales : num [1:1070] 34160 12438 47795 41250 37952 ...
## $ publisher revenue : num [1:1070] 20496 7462 28677 24750 22772 ...
## $ sale price : num [1:1070] 4.88 1.99 8.69 7.5 7.99 4.99 6.99 6.89 6.99 6.49 ...
## $ sales rank : num [1:1070] 1 2 3 3 4 5 6 8 9 10 ...
## $ Publisher : chr [1:1070] "HarperCollins Publishers" "HarperCollins Publishers" "Amazon Digital Services, Inc." "Hachette Book Group" ...
## $ units sold : num [1:1070] 7000 6250 5500 5500 4750 ...
## - attr(*, "spec")=
## .. cols(
## .. index = col_double(),
## .. `Publishing Year` = col_double(),
## .. `Book Name` = col_character(),
## .. Author = col_character(),
## .. language_code = col_character(),
## .. Author_Rating = col_character(),
## .. Book_average_rating = col_double(),
## .. Book_ratings_count = col_double(),
## .. genre = col_character(),
## .. `gross sales` = col_double(),
## .. `publisher revenue` = col_double(),
## .. `sale price` = col_double(),
## .. `sales rank` = col_double(),
## .. Publisher = col_character(),
## .. `units sold` = col_double()
## .. )
## - attr(*, "problems")=<externalptr>#convert author rating to factor
book_sales$Author_Rating <- factor(book_sales$Author_Rating)
# Check if Author_Rating is a factor variable
is.factor(book_sales$Author_Rating)## [1] TRUE#round integers
book_sales$Book_average_rating <- round(book_sales$Book_average_rating)
book_sales$sales_rank <- round(book_sales$`sales rank`)
book_sales$sale_price <- round(book_sales$`sale price`)
#generalized linear model with poisson distribution
glm_model <- glm(`gross sales` ~ Book_average_rating + `sales rank` + `sale price` + Author_Rating,
data = book_sales, family = poisson)
summary(glm_model)##
## Call:
## glm(formula = `gross sales` ~ Book_average_rating + `sales rank` +
## `sale price` + Author_Rating, family = poisson, data = book_sales)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.600e+00 3.346e-02 286.911 <2e-16 ***
## Book_average_rating -3.425e-01 8.314e-03 -41.195 <2e-16 ***
## `sales rank` -3.456e-03 3.376e-06 -1023.636 <2e-16 ***
## `sale price` 1.097e-01 1.138e-04 964.171 <2e-16 ***
## Author_RatingFamous -2.705e-03 4.730e-03 -0.572 0.567
## Author_RatingIntermediate 1.096e-01 2.342e-03 46.777 <2e-16 ***
## Author_RatingNovice 3.681e-01 9.125e-03 40.342 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3678126 on 1069 degrees of freedom
## Residual deviance: 1010821 on 1063 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 5#extract predicted values from the Poisson regression model
predicted <- predict(glm_model, type = "response")
#sscatterplot of observed vs. predicted gross sales
plot(book_sales$`gross sales`, predicted, pch = 20, col = "blue",
xlab = "Observed Gross Sales", ylab = "Predicted Gross Sales",
main = "Observed vs. Predicted Gross Sales")
#represents perfect predictions
abline(0, 1, col = "red")
legend("topleft", legend = "Perfect Prediction", col = "red", lty = 1, pch = NA)
Diagnosing
To diagnose, I have chose to do a Residuals vs. \(\hat{y}\) values plot.
# residuals vs. predicted values plot
plot(glm_model$fitted.values, residuals(glm_model),
xlab = "Predicted Values", ylab = "Residuals",
main = "Residuals vs. Predicted Values")
abline(h = 0, col = "red")
The plot does show widening of residuals as the predicted values increase, indicating the presence of heteroscedasticity. This shows a problem with my model, as utilizing a poisson regression model indicates constant variance, which is not currently present.