Exploring the response variable

Suppose we are considering the following model for price:

\[Y_i \stackrel{iid}{\sim} N(\mu,\sigma)\] This models assumes that the distribution of price is unimodal and symmetric.

When we find ourselves dealing with a skewed distribution, but we would like to use a normal model, we generally consider taking the natural log of the variable.

Based on what we seen explored so far, we will consider the model

\[log(Y_i) \stackrel{iid}{\sim} N(\mu,\sigma)\]

In this model, \(\mu\) represents the average log price of one night in an airbnb rental in Chicago. \(\sigma\) represents the standard deviation of log rental prices, or how much we expect the log price of an individual rental would differ from \(\mu\).

However, one of the “i”s in iid stands for independent. This suggests that the prices of each of the rentals can be treated as independent of one another. In other words, we essentially assume the rentals are a random sample from across Chicago and the same \(\mu\) and \(\sigma\) are a good description for all rentals.

If we look at our data set, however, we will notice something about the structure. The city of Chicago is divided into several districts, and district is a variable in a our data set. If we have multiple rentals per district in our data set, it is then not appropriate to consider each rental as independent.

knitr::kable( table(  ), col.names = c("District", "Number of Rentals"))

As we can see in the table, there are multiple rentals per district! This means that we have a grouping variable and we will incorporate this into the model using a random effect.

EDA

When we are considering adding random effects, it is helpful to explore how the variable we are modeling differs from group to group. To see how the log prices changes from district to district, we need to create a visualization that can look at the relationship between a categorical variable (district) and a numeric one (log price). A good plot for such a purpose is a side-by-side boxplot.

ggplot( bnb, aes( x= district, y = log(price))) + 
  geom_boxplot( )

These plots indicate that log prices vary from district to district. However, based on our model we will be modeling the average log price per district. Box plots show medians, not means.

To compute the district specific average log price, we want to (1) gather together all the rentals in each district and (2) find the mean log rental price for each district. One really helpful package for doing such multi-stage tasks in R is dplyr.

library(dplyr)

This package allows us to use a tool called a pipe (|>) which allows us to do multi-stage tasks. Take a look at the code below.

bnb |>
    group_by(district) |>
    summarise( mean( log(price) )) |>
    knitr::kable()

Hmm. This looks confusing. To make this easier to read, I am going to annotate the code. This means I will put in lines of text (not code) that will explain what each line of code is doing. To annotate in R, you just put a # in front of anything you want R to treat as text rather than code.

# Start with this data set AND THEN
bnb |>
  # Gather together all rentals in a district AND THEN 
  group_by(district) |>
  # Find the average log price for each district AND THEN
  summarise( mean(log(price) ) ) |> 
  # Format the output
  knitr::kable(col.names=c("District", "Avg. Log Price"))

Based on the annotated code above, we can see that the pipe |> can be read as “and then”. So, this code means start with the bnb data set (bnb) and then (|>) collect all the rentals in the same district together (group_by(district)) and then (|<) take the mean of the price in each district (summarise( mean(log(price) ) )).

Fitting the LME Model

Based on the work we have done so far, the first model we will consider for price is:

\[log(Price_{ij}) {\sim} N( \mu_j , \sigma_\epsilon),\]

\[\mu_j \stackrel{iid}{\sim} N( 0 , \sigma_u)\]

\[\epsilon_{ij} \stackrel{iid}{\sim} N( 0 , \sigma_\epsilon)\]

where \(i = 1, \dots, 1561\) denotes the rental and \(j = 1, \dots, 9\) denotes the district.

To fit this model in R, we need the lme4 library.

suppressMessages(library(lme4))

We can then use the following code::

LME1 <- lmer(  log(y) ~  (1|group), data = )

where

• y is the name of our y variable.
• group is the name of the variable corresponding to the random effect.
• data = should have the name of our data set after the =

Considering \(u_j\)

Suppose we want to see the different values of \(\hat{\beta}_0 + \hat{u}_j\) that we have estimated in our model. We have assumed that average log price differs from district to district, so we might be interested to see how much those district specific averages differ from the overall average.

To find that information, we use the code:

coef(LME1)$district

This output provides the log average price in each district. Note that district is our grouping variable; if you have a data set with a different grouping variable, you need to change district to match that new grouping variable.

Option 1 in R:

# The top part of the output
knitr::kable(summary(LME2)$varcor)

# The bottom part of the output
knitr::kable(summary(LME2)$coefficients)

To visualize the fitted model, we can use the following:

# Using x = bedrooms and y = log price
g1<- ggplot(bnb, aes(x = bedrooms, y = log(price), color = district)) +
  # Make a scatter plot 
  geom_point() +
  # Add labels 
  labs(y= "Bedrooms", x = "Log Price",
  title = "No Random Slope") +
  # Make the background of the plot clear
  theme_light() +
  # Remove the legend
  theme(legend.position = "none") +
  # Add in the lines 
  geom_abline(slope = .42 , intercept = coef(LME2)$district[,1])
  
# Show the plot 
g1

Option 2 in R

Right now, our model assumes that as the number of bedrooms increases by 1, the change in the log price is the same for every neighborhood. However, what if we think the change in log price if we increase the number of bedrooms might also change with district? In this case, we consider adding a random slope.

There is one part of the output from the LME model with a random slope that we have discussed yet. In the top part of the LME output, you will see a column called Corr, which stands for correlation. Specifically, this column tells us the correlation between our random intercept for district and the random slope for district.

What does this mean?? Well, a negative correlation means districts with a lower district-average tend to have steeper slopes than districts with higher district-averages. In other words, if you are in a district with lower rents, the number of bedrooms makes more of a difference in the rent.

The magnitude of hte correlation indicates the strength of this relationship. The closer to 1 or -1, the stronger the correlation.

We can visualize the output of this model by:

g2 <- ggplot(bnb, aes(x = bedrooms, y = log(price), color = district))+
     # Make a scatter plot 
     geom_point() +
     # Add labels 
     labs(y= "Bedrooms", x = "Log Price", title = "Random Slope") +
     # Make the background of the plot clear
     theme_light() +
     # Remove the legend
     theme(legend.position = "none") +
     # Add in the lines 
     geom_abline(slope = coef(LME3)$district[,2] , 
     intercept = coef(LME3)$district[,1])
     
# Show the plot 
g2 

Comparing the two options

At the moment, we have two different possibilities: a model with a random slope and a model with no random slope. Which should we choose??

We can look at plots of the two fitted models side by side using the following:

gridExtra::grid.arrange(g1,g2,ncol = 2)

If the plots look very different, this means that by including a random slope, we change the fitted model quite a bit. If the two plots look similar or only a little different, this means choosing to include a random slope results in a fitted model that is very similar to a model with no random slope.

In the case of little change, a random slope may not be needed. Random slopes drastically increase the difficulty in estimating the parameters in our LME model, so we don’t generally want to include them if they are not needed.

If we want to numerically compare the two options, we can use the AIC. This is another tool to help us decide whether to include the random slope.

AIC(LME2)
AIC(LME3)

Using the Final Model

Now that we have selected a final model, we can use it!

References