Introduction

Research Question: Does number of minifigures, number of pieces, number of instruction booklet pages, or year of release best predict the recommended retail price of different lego sets?

About the dataset:

We will be using the lego_population dataset from openintro.org. The data set has 1304 rows and 14 columns. The 6 columns relevant to our research question are: set_name, which has the retail name of each set; price, which is the MSRP of each set; pieces, which is the part count of each set; year, which is the year the set was first released; pages, which is the number of pages in the set’s instruction manual; and minifigures, which is the number of minifigures included in the set.

Data Analysis

First, let’s load in the data and start cleaning it by selecting only the columns relevant to us.

df <- read.csv('lego_population.csv')
cleaned_df <- df |>
    select(set_name, price, pieces, year, pages, minifigures)

Next, we will check that the data has loaded in properly. We are expecting 1304 rows and 6 columns. price should be numeric, and the other columns should be either numeric or integer type.

# check number of rows and columns
str(cleaned_df)
## 'data.frame':    1304 obs. of  6 variables:
##  $ set_name   : chr  "Extra Dots - Series 2" "Extra Dots - Series 1" "Creative Blue Bricks" "Creative Green Bricks" ...
##  $ price      : num  3.99 3.99 4.99 4.99 4.99 4.99 4.99 4.99 4.99 4.99 ...
##  $ pieces     : int  109 109 52 60 33 33 33 33 33 33 ...
##  $ year       : int  2020 2020 2020 2020 2020 2020 2020 2020 2020 2020 ...
##  $ pages      : int  NA NA 37 37 NA NA NA NA NA NA ...
##  $ minifigures: int  NA NA NA NA NA NA NA NA NA NA ...

Everything looks good. Next let’s look at the top and bottom of the dataset to ensure the data were read in properly, things are properly formatted, and that everything is there.

# check top and bottom of data set
head(cleaned_df)
##                   set_name price pieces year pages minifigures
## 1    Extra Dots - Series 2  3.99    109 2020    NA          NA
## 2    Extra Dots - Series 1  3.99    109 2020    NA          NA
## 3     Creative Blue Bricks  4.99     52 2020    37          NA
## 4    Creative Green Bricks  4.99     60 2020    37          NA
## 5   Funky Animals Bracelet  4.99     33 2020    NA          NA
## 6 Sparkly Unicorn Bracelet  4.99     33 2020    NA          NA
tail(cleaned_df)
##                           set_name  price pieces year pages minifigures
## 1299       Aurora's Forest Cottage  39.99    300 2020    NA          NA
## 1300               SPIKE Prime Set 329.95    528 2020    NA           2
## 1301         Mario's House & Yoshi  29.99    205 2020    NA           2
## 1302          Toad's Treasure Hunt  69.99    464 2020    NA           4
## 1303   Bowser's Castle Boss Battle  99.99   1010 2020    NA          NA
## 1304 Propeller Mario Power-Up Pack   9.99     13 2020    NA          NA

Everything looks okay, but we can already see many missing values. To make our multiple regression model, we need to either remove all rows with missing values or impute them. Let’s check for NA values in each column, and check how many total rows have missing data.

colSums(is.na(cleaned_df))
##    set_name       price      pieces        year       pages minifigures 
##           0         239          42           0         154         450
# Total rows with any missing data
sum(is.na(cleaned_df$price) | 
    is.na(cleaned_df$pieces) |
    is.na(cleaned_df$year) |
    is.na(cleaned_df$pages) |
    is.na(cleaned_df$minifigures))
## [1] 549
# 549 observations, or about 42% of the data

Removing 42% of the data would have an extreme effect on our results, so we will impute. We will make box plots of each column with missing values to see if they have any outliers.

#print # of outliers to see if it makes sense to impute with median

#price
length(boxplot.stats(cleaned_df$price)$out)
## [1] 88
#pieces
length(boxplot.stats(cleaned_df$pieces)$out)
## [1] 98
#pages
length(boxplot.stats(cleaned_df$pages)$out)
## [1] 60
#minifigs
length(boxplot.stats(cleaned_df$minifigures)$out)
## [1] 24

They all feature many outliers, so we will impute with median. It is worth noting that this much imputation will likely bias standard errors and shrink variance in our model.

cleaned_df <- cleaned_df |>
  mutate(price = if_else(is.na(price), median(price, na.rm = TRUE), price)) |>
  mutate(pieces = if_else(is.na(pieces), median(pieces, na.rm = TRUE), pieces)) |>
  mutate(pages = if_else(is.na(pages), median(pages, na.rm = TRUE), pages)) |>
  mutate(minifigures = if_else(is.na(minifigures), median(minifigures, na.rm = TRUE), minifigures))

#check if any NA remain
colSums(is.na(cleaned_df))
##    set_name       price      pieces        year       pages minifigures 
##           0           0           0           0           0           0

Next, let’s check how much this changed our summary statistics.

summary(cleaned_df)
##       set_name        price            pieces             year     
##  Length   :1304   Min.   :  1.99   Min.   :   1.00   Min.   :2018  
##  N.unique :1264   1st Qu.: 19.99   1st Qu.:  83.75   1st Qu.:2018  
##  N.blank  :   0   Median : 29.99   Median : 198.00   Median :2019  
##  Min.nchar:   2   Mean   : 43.35   Mean   : 383.28   Mean   :2019  
##  Max.nchar:  66   3rd Qu.: 49.99   3rd Qu.: 440.25   3rd Qu.:2020  
##                   Max.   :699.99   Max.   :6020.00   Max.   :2020  
##      pages          minifigures    
##  Min.   :   1.00   Min.   : 1.000  
##  1st Qu.:  36.00   1st Qu.: 2.000  
##  Median :  68.00   Median : 2.000  
##  Mean   :  90.04   Mean   : 2.703  
##  3rd Qu.: 112.00   3rd Qu.: 3.000  
##  Max.   :1527.00   Max.   :28.000
summary(select(df, set_name, price, pieces, year, pages, minifigures))
##       set_name        price            pieces            year     
##  Length   :1304   Min.   :  1.99   Min.   :   1.0   Min.   :2018  
##  N.unique :1264   1st Qu.: 14.99   1st Qu.:  79.0   1st Qu.:2018  
##  N.blank  :   0   Median : 29.99   Median : 198.0   Median :2019  
##  Min.nchar:   2   Mean   : 46.35   Mean   : 389.4   Mean   :2019  
##  Max.nchar:  66   3rd Qu.: 49.99   3rd Qu.: 455.0   3rd Qu.:2020  
##                   Max.   :699.99   Max.   :6020.0   Max.   :2020  
##                   NAs    :239      NAs    :42                     
##      pages          minifigures    
##  Min.   :   1.00   Min.   : 1.000  
##  1st Qu.:  32.00   1st Qu.: 1.000  
##  Median :  68.00   Median : 2.000  
##  Mean   :  92.99   Mean   : 3.074  
##  3rd Qu.: 124.00   3rd Qu.: 4.000  
##  Max.   :1527.00   Max.   :28.000  
##  NAs    :154       NAs    :450

And finally, let’s see if there are any noticeable trends year-by-year in the average price of a set:

cleaned_df |>
  group_by(year) |>
  summarise(avg_price = mean(price, na.rm = TRUE)) |>
  arrange(year)
## # A tibble: 3 × 2
##    year avg_price
##   <int>     <dbl>
## 1  2018      40.3
## 2  2019      44.9
## 3  2020      44.9

Mean LEGO set price does trend up with release year! We will see how significant this trend is when we make the model in our next step.

Regression Analysis

Now it is time to answer our research question of which variable best predicts LEGO set price. We will use a Multiple Linear Regression Model because our outcome variable price is continuous. The factors in our research question (number of pieces, number of minifigures, number of instruction booklet pages, and year of release) were selected because they can be reasonably expected to influence the price of a set, be that because of inflation and market timing or because they directly affect the amount of materials in the box.

model <- lm(price ~ pieces + year + pages + minifigures, cleaned_df)
summary(model)
## 
## Call:
## lm(formula = price ~ pieces + year + pages + minifigures, data = cleaned_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -134.88  -10.17   -3.76    7.27  325.63 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7.810e+02  1.704e+03  -0.458   0.6469    
## pieces       7.049e-02  1.765e-03  39.939  < 2e-16 ***
## year         3.910e-01  8.441e-01   0.463   0.6433    
## pages        6.051e-02  1.010e-02   5.989 2.73e-09 ***
## minifigures  8.703e-01  3.789e-01   2.297   0.0218 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.99 on 1299 degrees of freedom
## Multiple R-squared:  0.7834, Adjusted R-squared:  0.7827 
## F-statistic:  1174 on 4 and 1299 DF,  p-value: < 2.2e-16

Piece count appears to be the most reliable and statistically significant predictor of a LEGO set’s retail price.

Let’s look at the confidence intervals for our coefficients:

# confidence intervals
confint(model)
##                     2.5 %       97.5 %
## (Intercept) -4.124263e+03 2.562356e+03
## pieces       6.702664e-02 7.395146e-02
## year        -1.264864e+00 2.046922e+00
## pages        4.068487e-02 8.032595e-02
## minifigures  1.271125e-01 1.613582e+00

Model Assumptions and Diagnostics

We must check for 5 assumptions:

  1. Linearity
  2. Homoscedasticity
  3. Independence of observations
  4. Normality of residuals
  5. Multicollinearity

Homoscedasticity & Linearity

  • To check the assumption of Homoscedasticity, we will make a Residuals vs. Fitted plot. Ideally, this plot will have a random scatter of points centered along the horizontal zero line. This would confirm that the model has a linear relationship, variance between errors is equal, and there are no major outliers.
# resid v. fitted plot to check Linearity, Homoscedasticity
plot(model, which = 1)

The residuals vs. fitted plot shows an approximately linear relationship in most of the data, but reveals clear heteroscedasticity. The spread of the residuals gets larger at higher fitted values, creating a cone or funnel shape. This suggests that the homoscedasticity assumption has been violated, and the model is less accurate at predicting prices of higher-MSRP sets.

  • We can further confirm heteroscedasticity with a Scale-Location plot:
plot(model, which = 3)

The Scale-Location plot shows further evidence of heteroscedasticity. The red trend line goes higher as fitted values increase, showing that the model is indeed less accurate when predicting prices for high-priced sets.

Independence of observations

Because each sample is an individually released LEGO set with no obvious influence on the MSRP of other LEGO sets, it is reasonable to assume independence. However, the data contains many themed LEGO sets (Harry Potter, Star Wars, LEGO City, etc.) that may share pricing conventions with other sets in that theme. The model does not account for this.

Normality of residuals

  • To check the normality of residuals we will make a Q-Q plot. Ideally the points on the plot will create a straight diagonal line. An S-curve or heavy tails would indicate that the residuals are not normally distributed, which would make our p-values and confidence intervals less valid.
# Q-Q plot to check normality of residuals
plot(model, which = 2)

The Q-Q plot indicates that the normality assumption is violated. We can see that the Q-Q plot has a rather extreme tail on the right-hand side, indicating that the residuals skew right rather than having a normal distribution. The right-hand tail has three extreme outliers (observations 826, 1300, 1240).

No Multicollinearity

  • To check for multicollinearity, we will make a correlation matrix and check for highly correlated predictors
library(corrplot)
## corrplot 0.95 loaded
# select only predictors
cor_matrix <- cor(select(cleaned_df, pieces, year, pages, minifigures))
corrplot(cor_matrix, method = "number")

pieces and pages have a fairly strong correlation of 0.73, indicating multicollinearity. This means that the standard errors of those two predictors’ coefficients may be inflated. minifigures and pieces; as well as minifigures and pages have some correlation with scores of 0.46 and 0.41 respectively. These scores are less dramatic but still worth noting.

Residuals vs. Leverage Plot

  • We know there are some extreme outlier sets in our data: observations 826, 1300, and 1240. They appear labeled in our Q-Q, Scale-Location, and Residuals vs. Fitted plots. Our model had an unusually difficult time predicting the retail prices of these sets. Let’s find their names and see if there is anything special about them.
  cleaned_df |>
  filter(row_number() %in% c(826, 1300, 1240))
##                  set_name  price pieces year pages minifigures
## 1 Imperial Star Destroyer 699.99   4784 2019   444           2
## 2             My XL World 239.95    480 2020    68           2
## 3         SPIKE Prime Set 329.95    528 2020    68           2

Finding each set online reveals that each one is a unique oddity.

  • The Imperial Star Destroyer is a licensed set with an intricate internal design that yields a very physically large toy ship. These three cost-increasing factors are not accounted for in the model, so it would have no way of predicting this set’s cost accurately.

  • The SPIKE Prime Set is a robotics kit sold directly to secondary schools through a program called “LEGO Education.” It includes motors and sensors that are much more costly to produce than an average LEGO brick, as well as components that do not contribute to part count such as sorting trays and computer programming software.

  • The My XL World set is another LEGO Education set for children ages 2-6. The set has oversized bricks that use more plastic than typical smaller LEGO parts, and our model has no way of telling the difference.

These sets have extreme residuals. That means that their actual price is unusually far from the price predicted by the model. To see how they affect the fit of the model, we can make a Residuals vs Leverage plot and see the influential cases (sets that, if omitted, would change the fit of the model):

# resid. vs leverage plot
plot(model, which = 5)

* The residuals vs. leverage plot shows that most sets have low leverage, indicating expected predicted prices from the model. Points 815 and 822 have notably higher leverage, which could mean they have very high piece or page counts, but their residuals are fairly small, meaning the model still predicts their prices reasonably well. Point 826, the Imperial Star Destroyer, has an very high residual (~14) and a leverage higher than most sets (~0.06). It is within the Cook’s distance 0.5 contour on the graph but sits close to the border, meaning that it has a moderate influence on the model without crossing the conventional threshold for an influential case.

Conclusion & Future Directions

Between number of minifigures included, number of pieces in the set, number of pages in the instruction booklet, and year of release, the number of pieces in the set is the strongest, most significant predictor of LEGO set price. pages and minifigures are also significant, but with smaller effects. year is not a significant predictor when the other three variables are controlled for. Our model explains a substantial ~78% of variance in LEGO set price with a Multiple \(R^2\) of 0.7834 and an Adjusted \(R^2\) score of 0.7827.

We can conclude that piece count best predicts LEGO set price based on its statistical significance compared to the other predictors. This suggests that LEGO’s pricing is primarily driven by the physical quantity of bricks in the set.

Limitations

Based on our diagnostics, we can observe many limitations in our model:

  • Heteroscedasticity: Our model is less reliable for expensive sets

  • Non-normal residuals: Our residuals are very right-skewed rather than strictly normally distributed. This is largely because of our three extreme outliers: the Star Wars ship and the two LEGO Education sets.

  • Multicollinearity: pieces and pages have a strong correlation (r = 0.73). This makes their coefficients in the model less reliable.

  • 42% of rows contain imputed medians: This is likely reducing variance and biasing standard error since so many values are the same. The minifigure standard error likely got the worst of this, as that column saw the most imputations.

Future Directions

Because of the limitations in my knowledge and the number of columns that went unused in this project, there is much more to explore with this dataset.

  • One could refit the model excluding the Imperial Star Destroyer (fairly influential according to Residual-Leverage plot) to see how it changes

  • One could add a theme variable as a predictor or a grouping factor to see how the relationship between piece and price for licensed sets (Star Wars, Marvel) compares to that of original LEGO IPs.

  • One could make a new model without the page predictor to address collinearity concerns.

  • The education set with the large bricks had an extremely high residual. One could fit a multiple regression using brick size and part count as factors to predict price. This article explains an activity like this for students learning linear regression.

References:

https://stackoverflow.com/questions/29787850/how-do-i-add-a-url-to-r-markdown

https://www.tandfonline.com/doi/full/10.1080/26939169.2021.1946450

https://www.geeksforgeeks.org/machine-learning/r-squared-vs-adjusted-r-squared-difference/ (Multiple \(R^2\) vs. Adj. \(R^2\))

https://library.virginia.edu/data/articles/diagnostic-plots (interpreting Residuals vs. Leverage plots)

https://statsandr.com/blog/outliers-detection-in-r/ (Detecting outliers)

LEGO sets:

https://brickset.com/sets/45678-1/SPIKE-Prime-Set

https://brickset.com/sets/45028-1/My-XL-World

https://brickset.com/article/45525/review-75252-imperial-star-destroyer