Linear Modeling in R

Quick Intro to R Notebooks and R Markdown

This is an R Markdown Notebook. When you execute R code within the notebook, the results appear beneath the code.

This file was created in RStudio by going to File…New File…R Notebook.

R code needs to be in “chunks” in an R Markdown Notebook. Below is an example of an R code chunk. It makes a parabola.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter (Win/Linux) or Cmd+Shift+Return (Mac).

x <- seq(-1, 1, by = 0.01)
y <- x^2
plot(x, y, type = "l")

To hide the output, click the Expand/Collapse output button. To clear results (or an error), click the “x”.

You can also press Ctrl+Enter (Win/Linux) or Cmd+Return (Mac) to run one line of code at a time (instead of the entire chunk).

Add a new R code chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I (Win/Linux) or Cmd+Option+I (Mac).

CODE ALONG 0

Insert a new R code chunk below and type and run the code: Sys.time()

Linear Modeling with Simulated Data

Instead of using theory and formulas, let’s explore and review linear modeling using simulated data.

Below we assign to x the values 1 - 25. Then we generate y as a function of x using the formula 10 + 5*x:

x <- 1:25
y <- 10 + 5*x  # formula for a line
d <- data.frame(x, y)
plot(y ~ x, data = d)

• 10 is the intercept.
• 5 is the slope.
• y is completely determined by x (y = 10 + 5*x)

Now let’s add some “noise” to our data by adding random draws from a Normal distribution with mean = 0 and a standard deviation = 10. The rnorm() function allows us to draw random values from a Normal distribution.

set.seed(1) ensures we all generate the same “random” data:

set.seed(1)
noise <- rnorm(n = 25, mean = 0, sd = 10)
# Add the noise to 10 + 5*x and re-draw plot
d$y <- 10 + 5*x + noise
plot(y ~ x, data = d)

Now y appears to be associated with x, but not completely determined by x.

y is the combination of a fixed part and a random part:

1. fixed: 10 + 5*x
2. random: rnorm(n = 25, mean = 0, sd = 10)

What if we were given this data and told to determine the process that generated it? In other words, work backwards and fill in the blanks:

1. ___ + ___*x
2. rnorm(n = 25, mean = 0, sd = ____)

That’s one way to think of linear modeling/regression. You have some numeric response (or dependent) variable, and you want to find the model (or the formula) that generated the data.

Traditional linear modeling/multiple regression assumes the following (among others):

1. the formula is a weighted sum of predictors (eg, y = 10 + 5*x)
2. the noise is a random draw from a Normal distribution with mean = 0
3. the standard deviation of the Normal distribution is constant (eg, 10)

Linear modeling tries to recover the weights in the first assumption and the standard deviation in the third assumption.

Let’s demonstrate. Below we attempt to recover the data generating process for our data. For this we use the lm() function. We have to specify the formula for the first assumption. The 2nd and 3rd assumptions are built into lm().

The formula “y ~ x” means we think the model is “y = intercept + slope*x”. (Unless we specify otherwise, this assumes we want to estimate an intercept.) This tells lm() to take our data and find the best intercept and slope. Notice this is the correct model!

When you use lm() you’ll usually want to save the results to an object. Below I save it to “mod”. Then we view the results of the model using summary()

mod <- lm(y ~ x, data = d)
summary(mod)

Call:
lm(formula = y ~ x, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-23.876  -4.613   2.254   5.909  14.648 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   11.135      3.999   2.784   0.0105 *  
x              5.042      0.269  18.743 1.98e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.7 on 23 degrees of freedom
Multiple R-squared:  0.9386,    Adjusted R-squared:  0.9359 
F-statistic: 351.3 on 1 and 23 DF,  p-value: 1.981e-15

The model returns the following estimates:

1. y = 11.135 + 5.042 * x
2. noise = rnorm(n = 25, mean = 0, sd = 9.7)

These are pretty close to the “true” values of 10, 5, and 10 we used to generate the data.

In real life, we DO NOT KNOW the formula in part 1. The real data generation process will be far more complicated. The formula we propose will just be an approximation and may not be good.

In real life, we DO NOT KNOW if the Normality assumption or constant variance assumption of the noise is plausible.

How can we evaluate our model formula?

We could use our model estimates to generate data and see if they look similar to our original data. Run entire chunk at once and run more than once. The black points don’t change but the red ones do. That looks pretty good! Our model-generated data appears similar to our observed data.

# d$y <- 10 + 5*x + noise
d$y2 <- 11.135 + 5.042*d$x + rnorm(25, 0, 9.7)
plot(y ~ x, data = d) # original data
points(d$x, d$y2, col = "red") # simulated data

We can also compare smooth density curves of the original and model-generated data. Smooth density curves are basically smooth versions of histograms. If we have a good model, data generated by our model should have a similar distribution to the original data. Run entire chunk at once and run more than once.

hist(d$y, freq = FALSE) # freq = FALSE means area of bars sums to 1
lines(density(d$y))  # original data

d$y2 <- 11.135 + 5.042*d$x + rnorm(25, 0, 9.7)
lines(density(d$y2), col = "red")  # simulate data

This looks good as well. The distribution of our model-generated data is very similar to our observed data. You should do this more than once, say 50 times, to ensure the model consistently generates data similar to the observed data. We show one way to do that later in the workshop.

Since we think our model is “good”, we might use it to make a prediction. For example, when x = 10 what’s the expected value of y? Put another way, what’s the mean of y conditional on x = 10? We can do that with the predict() function. The interval = "confidence" arguments says return a 95% confidence interval (CI) for this mean.

predict(mod, newdata = data.frame(x = 10), interval = "confidence")

       fit     lwr      upr
1 61.55932 57.2126 65.90603

The expected mean of y when x = 10 is about 61.6 with a 95% CI of (57.2, 65.9). The CI gives us some notion of how uncertain this expected mean is. In fact it would be better to report this as, “the expected mean of y when x = 10 is about 57 to 66.”

We might also try to summarize the relationship between y and x by examining the coefficients (or weights) in the summary output. We can extract the coefficients from the summary with the coef() function:

coef(summary(mod))

             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 11.134859  3.9994866  2.784072 1.054874e-02
x            5.042446  0.2690346 18.742742 1.981382e-15

The x coefficient says that y increases by about 5 for every one-unit increase in x, give or take about 0.27. The standard error gives us some indication of the uncertainty in this estimate. We talk more about t values and p values below.

TO SUMMARIZE: This is basic linear modeling:

1. propose and fit a model
2. determine if the model is good and that assumptions are mostly met
3. use the model to explain relationships and/or make predictions

CODE ALONG 1

Let’s see what happens when we fit a bad model. Below we add a new column to data frame d named z, which is a random sample of numbers on the range of -100 to 100. runif() samples numbers from a uniform distribution.

set.seed(4)
d$z <- runif(25, min = -100, max = 100)

REMINDER: Add a new R code chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I (Win/Linux) or Cmd+Option+I (Mac).

1. Model y as a function of z using lm(y ~ z, data = d) and save to an object called mod2. View the summary. What formula does the model return? What is the estimated standard deviation of the Normally distributed noise?

2. Use the model to simulate density histograms and compare to the original density histogram of d$y. Run the chunk several times to see how the model-generated density curve varies.

Let’s do some linear modeling with real data.

Import data

Let’s import the data we’ll be using today. The data we’ll work with is Albemarle County real estate data which was downloaded from the Office of Geographic Data Services. We’ll use a random sample of the data.

URL <- 'https://raw.githubusercontent.com/clayford/dataviz_with_ggplot2/master/alb_homes.csv'
homes <- read.csv(file = URL)
homes$hsdistrict <- factor(homes$hsdistrict)
homes$cooling <- factor(homes$cooling)

Let’s look at the first few rows:

head(homes)

Variable name definitions:

• yearbuilt: year house was built
• finsqft: size of house in number square feet
• cooling: ‘Central Air’ versus ‘No Central Air’
• bedroom: number of bedrooms
• fullbath: number of full bathrooms (toilet, sink and bath)
• halfbath: number of half bathrooms (toilet and sink only)
• lotsize: size of land on which home is located, in acres
• totalvalue: total assessed value of home and property
• esdistrict: the elementary school the home feeds into
• msdistrict: the middle school the home feeds into
• hsdistrict: the high school the home feeds into
• censustract: the census tract the home is located in
• age: of the house in years as of 2018
• condition: assessed condition of home (Substandard, Poor, Fair, Average, Good, Excellent)
• fp: indicator if house has fireplace (0=no, 1=yes)

Linear Modeling with Real Estate Data

Let’s say we want to model the mean total value of a home as a function of various characteristics such as lot size, finished square feet, presence of central air, etc.

Let’s see how total value is distributed using a histogram. Notice it’s very skewed.

hist(homes$totalvalue)

We can also create a smooth density plot, which is a smooth version of a histogram:

plot(density(homes$totalvalue))

In order to model mean totalvalue of homes as a function of various characteristics, we need to propose a linear model. Unlike the previous example this is not simulated data for which we know the data generating process. How to propose a model? It helps to have some subject matter expertise.

Let’s fit a linear model using finsqft, bedrooms and lotsize. The plus (+) sign means “include” in model.

m1 <- lm(totalvalue ~ finsqft + bedroom + lotsize, data = homes)
summary(m1)

Call:
lm(formula = totalvalue ~ finsqft + bedroom + lotsize, data = homes)

Residuals:
     Min       1Q   Median       3Q      Max 
-1164152   -82296    -7690    57164  5879188 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -133328.25   16273.88  -8.193 3.73e-16 ***
finsqft         284.46       5.37  52.967  < 2e-16 ***
bedroom      -13218.41    5750.70  -2.299   0.0216 *  
lotsize        4268.77     219.32  19.464  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 227200 on 3021 degrees of freedom
Multiple R-squared:  0.6164,    Adjusted R-squared:  0.616 
F-statistic:  1618 on 3 and 3021 DF,  p-value: < 2.2e-16

The coef() function extracts the coefficients (or weights):

coef(m1)

 (Intercept)      finsqft      bedroom      lotsize 
-133328.2482     284.4613  -13218.4091    4268.7655 

This translates to:

totalprice = -133328.2482 + 284.4613*finsqft + -13218.4091*bedroom +
               4268.7655*lotsize`

Some naive interpretation:

• each additional finished square foot adds about $284 to price
• each additional bedroom drops the price by $13,218 (?)
• each additional lot size acre adds about $4268 to price

Each of these interpretations assumes all other variables are held constant! So adding a bedroom to a house, without increasing the lot size or finished square feet of the house, is estimated to drop the value of a home. Does this make sense?

Is this a “good” model? Let’s simulate data from the model and compare it to our observed data. A “good” model should generate data that looks similar to the original data.

We could do this by hand:

sim_values <- -133328.2482 + 284.4613*homes$finsqft + 
  -13218.4091*homes$bedroom + 4268.7655*homes$lotsize + 
  rnorm(3025, mean = 0, sd = 227200)

An easier and faster way is to use the simulate() function which allows you to generate multiple samples. Here we generate 50 samples. Each sample will have the same number of observations as our original sample (n = 3025). Each sample value is generated using our observed values for finsqft, bedroom, and lotsize. The result is a data frame with 50 columns.

sim1 <- simulate(m1, nsim = 50)

Now let’s plot our simulated data with our observed data using smooth density plots. We use a for loop to add smooth density estimates of the 50 simulations. The syntax sim1[[i]] extracts column i as a vector. (Run entire chunk at once.)

plot(density(homes$totalvalue))
for(i in 1:50)lines(density(sim1[[i]]), col = "grey80")

(See end of notebook for how to create this plot using ggplot2 and for how to turn this into a function.)

This does not appear to be a good model. In fact some of our simulated values are negative!

Before we revise our model recall the main assumptions:

1. totalvalue can be modeled by a weighted sum: totalvalue = Intercept + finsqft + bedrooms + lotsize
2. noise/error is from Normal dist’n with mean 0
3. the SD of the Normal dist’n is constant

R provides some basic diagnostic plots to assess 2 and 3. Just call plot on your model object

plot(m1)

How to interpret plots:

1. Residuals vs Fitted: for checking constant variance; should have a horizontal line with uniform and symmetric scatter of points; if not, evidence that SD is not constant

2. Normal Q-Q: for checking normality of residuals; points should lie close to diagonal line; if not, evidence that noise is not drawn from N(0, SD). This assumption is only really critical if you’re planning to use your model to predict exact values (instead of means) and calculate a confidence interval.

3. Scale-Location: for checking constant variance; should have a horizontal line with uniform scatter of point; (similar to #1 but easier to detect trend in dispersion)

4. Residuals vs Leverage: for checking for influential observations; points outside the contour lines are influential observations. Leverage is the distance from the center of all predictors. An observation with high leverage has substantial influence on the fitted value.

By default the 3 “most extreme” points are labeled by row number. 2658 appears in all four plots. It’s a really big expensive home.

homes[2658,]

These plots reveal that our assumptions of normally distributed residuals and constant variance are highly suspect. Our model is just bad.

What can we do?

Non-constant variance can be evidence of a wrong model or a very skewed response (or a bit of both). Recall that our response is quite skewed:

hist(homes$totalvalue)

When dealing with a response that is strictly positive and very skew (like dollar amounts), it is common to transform the response to a different scale. A common transformation is a log transformation. When we log transform totalvalue, the distribution looks a little more symmetric, though it’s important to note that is not an assumption of linear modeling!

hist(log(homes$totalvalue))

Let’s try modeling log-transformed totalvalue.

m2 <- lm(log(totalvalue) ~ finsqft + bedroom + lotsize, data = homes)

The diagnostic plots look better.

plot(m2)

But is this a “good model”? Is our proposed model of weighted sums good? Let’s simulate data and compare to observed data.

sim2 <- simulate(m2, nsim = 50)
plot(density(log(homes$totalvalue)))
for(i in 1:50)lines(density(sim2[[i]]), lty = 2, col = "grey80")

This doesn’t look too bad!

Let’s say we’re happy with this model. How do we interpret the coefficients? Since the response is log transformed, we interpret the coefficients as multiplicative instead of additive. Below we view the coefficients rounded to 4 decimal places.

round(coef(m2), 4)

These are proportions. To get percentages, multiply by 100.

round(coef(m2), 4) * 100

Some naive interpretations:

• each additional finished square foot increases expected mean price by 0.05%. Or multiply by 100 to say each additional 100 finished square feet increases expected mean price by 5%.
• each additional bedroom increases expected mean price by about 4.3%
• each additional lot size acre increases expected mean price by about 0.47%

Remember, the interpretation assumes all other variables are held constant!

A slightly more precise estimation can be obtained by exponentiating the coefficients. Below we exponentiate using the exp() function and then round to 4 decimal places.

round(exp(coef(m2)), 4) 

For example, each additional bedroom (assuming all else held constant) increases expected total price by about 4.4%. Multiplying by 1.0439 is equivalent to adding 4.39%.

Let’s review the summary output:

summary(m2)

SUMMARY OVERVIEW

• Residuals section: quick assessment of residuals. Ideally 1Q/3Q and Min/Max will be roughly equivalent in absolute value.
• Coefficients: lists the estimated coefficients along with hypothesis tests for the null that each coefficient is 0. Est/SE = t-value.
• Residual standard error: estimate of the constant standard deviation of the normal dist’n of the residuals (noise)
• degrees of freedom: sample size - number of coefficients (3025 - 4)
• R-squared: proportion of variance explained
• F-statistic: overall test that all coefficients (except intercept) are 0.

All of the p-values refer to hypothesis tests that the coefficients are 0. Many statisticians and researchers prefer to look at confidence intervals.

round(confint(m2) * 100, 4)

According to our model, each additional bedroom adds about 3% to 6% to the value of a home, assuming all else held constant.

CODE ALONG 2

1. Insert a code chunk below and model log(totalvalue) as function of fullbath and finsqft. Call your model m3

2. Insert a code chunk below and check the diagnostic plots

3. How do we interpret the fullbath coefficient?

4. Insert a code chunk below and simulate data from the model and compare to the observed totalvalue. Does this look like a good model?

Categorical predictors

Let’s add hsdistrict to the model we just fit. Does being in a particular high school district affect total value of a home? Here’s a quick tally:

table(homes$hsdistrict)

These levels are not numbers, so how does R handle this in a linear model? It creates a contrast. This is a matrix of zeroes and ones. If you have K levels, you’ll have K-1 columns. In this case we’ll have two columns: one for Monticello HS and one for Western Albemarle HS. By default R takes whatever level comes first alphabetically and makes it the baseline or reference level.

Let’s look at the default contrast, called a treatment contrast. To do this we convert the hsdistrict column to factor and then use the contrasts() function. (We don’t have to do this to add hsdistrict to our model! I’m just doing this to output the contrast “definition”)

homes$hsdistrict <- factor(homes$hsdistrict)
contrasts(homes$hsdistrict)

A model with hsdistrict will have two coefficients: Monticello and Western Albemarle

• a home in Albemarle HS district gets two zeroes
• a home in Monticello HS district gets a one in the Monticello column
• a home in Western Albemarle HS district gets a one in the West Alb column

Let’s fit our new model.

m4 <- lm(log(totalvalue) ~ fullbath + finsqft + hsdistrict, data = homes)
summary(m4)

The coefficients for Monticello and Western Albemarle are in relation to Albemarle HS.

round(coef(m4) * 100, 4)

It appears that the value of a home in Western Albemarle will be about 10% higher than an equivalent home in Albemarle. Likewise it appears that the value of a home in the Monticello district will be about 7% less than an equivalent home in the Albemarle district.

CODE ALONG 3

1. Insert a code chunk below and model log(totalvalue) as function of fullbath, finsqft and cooling. Call your model m5.

2. What is the interpretation of cooling?

Modeling Interactions

In our model above that included hsdistrict we assumed the effects were additive. For example, it didn’t matter what high school district your home was in, the effect of fullbath or finsqft was the same. It also assumed the effect of each additional fullbath was the same regardless of how big the house was, and vice versa. This may be too simple.

Interactions allow the effects of variables to depend on other variables. Again subject matter knowledge helps with the proposal of interactions. As we’ll see interactions make your model more flexible but harder to understand.

R makes it simple to include interactions in models. Just indicate an interaction between two variables by placing a colon (:) between them. Below we include 2-way interactions. (You can have 3-way and higher interactions but they’re very difficult to interpret.)

m6 <- lm(log(totalvalue) ~ fullbath + finsqft + hsdistrict + 
           fullbath:finsqft + fullbath:hsdistrict + 
           finsqft:hsdistrict, data = homes)
summary(m6)

Interpretation is now much more difficult. We cannot directly interpret the main effects of fullbath, finsqft or hsdistrict. They interact. What’s the effect of finsqft? It depends on fullbath and hsdistrict.

Are the interactions “significant” or necessary? We can use the anova() function to evaluate this question. This runs a series of partial F-tests. Each row below is a hypothesis test. The null is the model with this predictor is the same as the model without the predictor. The anova tests below use what are called Type I sums of squares. This respects the order of the variables in the model. So…

• the first test compares a model with just an intercept to a model with an intercept and fullbath.
• the second test compares a model with an intercept and fullbath to a model with an intercept, fullbath and finsqft.
• And so on.

If the null is true, the F value should be close to 1.

anova(m6)

The interaction finsqft:hsdistrict doesn’t appear to contribute much to the model given that everything else above it is in the model.

Just because an interaction is significant doesn’t necessarily mean it’s interesting or worthwhile. Nor can we infer anything about the nature of the interaction from the ANOVA table.

Effect plots can help us visualize and make sense of models with interactions. Let’s make one using the ggeffects package and talk about what it’s showing.

library(ggeffects)
plot(ggpredict(m6, terms = c("fullbath", "hsdistrict")))
# place fullbath on x-axis, group by hsdistrict

What’s the effect of fullbath? It depends. It’s more dramatic in Western Albemarle and Monticello. Of course a lot of the difference comes at extreme values of fullbath. The “ribbons” around the lines are 95% confidence intervals.

What exactly was plotted? We can see by calling ggpredict without plot

ggpredict(m6, terms = c("fullbath", "hsdistrict"))

ggpredict used our model to make totalvalue predictions for various values of fullbath in the three school districts, holding finsqft at 1828 (the median of finsqft).

We can specify our values if we like. For example, make an effect plot for 1 - 5 bathrooms and hold finsqft at 2000:

plot(ggpredict(m6, terms = c("fullbath[1:5]", "hsdistrict"), 
          condition = c(finsqft = 2000)))

What about the effects of finsqft and fullbath? This is an interaction of two numeric variables. The second variable has to serve as a grouping variable when creating an effect plot. Below we set fullbath to take values 2 - 5 and finsqft to take values of 1000 - 4000 in steps of 500.

plot(ggpredict(m6, terms = c("finsqft[1000:4000 by=500]", "fullbath[2:5]")))

The effect of finsqft seems to taper off the more fullbaths a house has. But there are few large homes with 2 full baths, and likewise, few small homes with 5 full baths. Even though the interaction is “significant” in the model, it’s clearly a very small interaction and probably not important.

CODE ALONG 4

1. Insert a code chunk below and model log(totalvalue) as function of fullbath, finsqft, cooling, and the interaction of finsqft and cooling. Call your model m7. Is the interaction warranted?

2. Visualize the interaction using the ggpredict function. Perhaps use [1000:4000 by=500] to set the range of finsqft on the x-axis. How notable is this interaction?

Nonlinear Effects

So far we have assumed that the relationship between a predictor and the response is linear (eg, for a 1-unit change in a predictor, the response changes by a fixed amount). That assumption can sometimes be too simple and not realistic. Fortunately there are ways to fit non-linear effects in a linear model.

Here’s a quick example of simulated non-linear data: a 2nd degree polynomial.

x <- seq(from = -10, to = 10, length.out = 100)
set.seed(3)
y <- 1.2 + 2*x + 0.9*x^2 + rnorm(100, mean = 0, sd = 10)
nl_dat <- data.frame(y, x)
plot(y ~ x, nl_dat)

Clearly a straight-line model will not work well for this data. The relationship between x and y is not adequately summarized by a straight line model.

If we wanted to try to “recover” the weights we used in simulating this data we could fit a polynomial model using the poly() function in the formula syntax:

## EXAMPLE CODE; NOT INTENED TO BE RUN
nlm1 <- lm(y ~ poly(x, degree = 2, raw = TRUE), data = nl_dat)

However, the recommended approach to fitting non-linear effects is to use natural splines instead of polynomials. Natural splines essentially allow us to fit a series of cubic polynomials connected at knots located in the range of our data.

The easiest option is to use the ns() function from the splines package, which comes installed with R. ns stands for “natural splines”. The second argument is the degrees of freedom (df). It may help to think of df as the number of times the smooth line changes directions.

Frank Harrell states in his book Regression Model Strategies that 3 to 5 df is almost always sufficient. His basic advice is to allocate more df to variables you think are more important.

Let’s see how it works with our simulated data.

library(splines)
nlm2 <- lm(y ~ ns(x, df = 2), data = nl_dat)
summary(nlm2)

The summary output is impossible to interpret. Let’s visualize the fit with an effect plot.

library(ggeffects)
plot(ggpredict(nlm2, terms = "x"), show_data = TRUE)

Let’s return to the homes data and fit a non-linear effect for finsqft using a natural spline with 5 df. Below we also include hsdistrict and lotsize and allow finsqft and hsdistrict to interact.

nlm3 <- lm(log(totalvalue) ~ ns(finsqft, 5) + lotsize + hsdistrict +
           ns(finsqft, 5):hsdistrict, 
         data = homes)

The anova function allows us to assess the non-linear effect and the interaction. Some sort of interaction between finsqft and hsdistrict appears to be present.

anova(nlm3)

The summary output is impossible to interpret.

summary(nlm3)

Effect plots are our only hope for understanding this model. Below we plot predicted totalvalue over finsqft ranging from 1000 - 3000, grouped by hsdistrict.

plot(ggpredict(nlm3, terms = c("finsqft[1000:3000 by=250]", "hsdistrict")))

The effect of finsqft on totalvalue seems more dramatic in Western Albemarle once you go beyond 1500 sq ft.

Does the model simulate data similar in distribution to the observed data?

sim4 <- simulate(nlm3, nsim = 50)
plot(density(log(homes$totalvalue)))
for(i in 1:50)lines(density(sim4[[i]]), col = "grey80")

We should still check model assumptions.

plot(nlm3)

Homes 12, 40, 963 and 1810 seem to stand out. Let’s have a look.

h <- c(12, 40, 963, 1810)
homes[h,c("totalvalue", "finsqft", "lotsize")]

Homes 12 and 40 are very low in totalvalue and the model way overpredicts their values. Home 963 has a massive totalvalue with 0 acres of lotsize. Home 1810 is on 611 acres and that value has a high leverage on its fitted value.

cbind(observed = homes$totalvalue[h], fitted = exp(fitted(nlm3)[h]))

CODE ALONG 5

1. Insert a code chunk below and model log(totalvalue) as function of finsqft with a natural spline of 5 df, cooling, and the interaction of cooling and finsqft (natural spline of 5 df). Call your model nlm4.

2. Use the anova function to check whether the interaction appears necessary. What do you think?

3. Create an effect plot of finsqft by cooling. Maybe try [1000:5000 by=250] for the range of values for finsqft.

Wrap-up

This was meant to show you the basics of linear modeling in R. Hopefully you have a better grasp of how linear modeling works. Scroll down for a few more topics.

What we did today works for independent, numeric outcomes. We had one observation per home and our response was totalvalue, a number. Our models returned expected mean total value given various predictors. This is a pretty simple design.

Things get more complicated when you have, say, binary responses or multiple measures on the same subject (or home). A non-exhaustive list of other types of statistical modeling include:

• generalized linear models (for binary and count responses)
• multinomial logit models (for categorical responses)
• ordered logit models (for ordered categorical responses)
• mixed-effect or longitudinal linear models (for responses with multiple measures on the same subject or clusters of related measures)
• survival models (for responses that measure time to an event)
• time series models (for responses that exhibit, say, seasonal variation over time)

References

• Faraway, J. (2005). Linear Models in R. London: Chapman & Hall.
• Fox, J. (2002). A R and S-Plus Companion to Applied Regression. London: Sage.
• Harrell, F. (2015). Regression Modeling Strategies (2nd ed.). New York: Springer.
• Kutner, M., et al. (2005). Applied Linear Statistical Models (5th ed.). New York: McGraw-Hill.
• Maindonald J., Braun, J.W. (2010). Data Analysis and Graphics Using R (3rd ed.). Cambridge: Cambridge Univ Press.

Bonus material/topics cut for time

How does lm() work?

lm() estimates regression coefficients using ordinary least squares, or OLS. The formula for this using matrix algebra is expressed as follows:

\[\hat{\beta} = (X'X)^{-1}X'Y \]

where X is the model matrix (our predictors as they appear in the model) and Y is the dependent variable. The prime symbol ' means transpose the matrix. The -1 means take the inverse. R was developed to make calculations such as these quite easy.

Let’s revisit model m2.

formula(m2)

Here are the coefficients returned by lm():

round(coef(m2), 4)

Let’s find these using the formula above. We can use model.matrix() to get X, t() to transpose, and solve() to take the inverse. Perform matrix multiplication with %*%

X <- model.matrix(~ finsqft + bedroom + lotsize, data = homes)
Y <- log(homes$totalvalue)
B <- solve(t(X) %*% X) %*% t(X) %*% Y
round(B, 4)

While this is theoretically what lm() does, it actually uses more sophisticated methods for faster performance and protection against numeric instability. This page goes into more details if you’re interested in learning more:

https://genomicsclass.github.io/book/pages/qr_and_regression.html

ANOVA revisited

We can also use the anova() function to compare nested models. This means one model is a subset of another. For example, below we fit progressively more complicated models, building up to our model, m2: log(totalvalue) ~ finsqft + bedroom + lotsize

mod_00 <- lm(log(totalvalue) ~ 1, data = homes) # intercept-only
mod_01 <- lm(log(totalvalue) ~ finsqft, data = homes)
mod_02 <- lm(log(totalvalue) ~ finsqft + bedroom, data = homes)
mod_03 <- lm(log(totalvalue) ~ finsqft + bedroom + lotsize, data = homes)

The anova() function allows us to compare these nested models. The null hypothesis is that two models are the same, in the sense they explain the same amount of variance in the response variable, log(totalvalue). A low p-value provides evidence that the more complicated model with more variables is a better model. The usual way to use anova() to compare models is to list the smaller models first. Below are three tests:

1. Model 2 vs Model 1
2. Model 3 vs Model 2
3. Model 4 vs Model 3

The end result is that model 4 is superior to the other three models.

anova(mod_00, mod_01, mod_02, mod_03)

Notice this identical to calling anova() on the full model.

anova(m2)

Another approach is using Type II sums of squares, where each variable is tested assuming all other variables are in the model. One approach to performing this test is the drop1() function. The name comes from the fact we’re dropping one variable at a time. The null hypothesis is dropping the variable from the model has no effect. A low p-value provides evidence against this hypothesis. Below each are three tests:

1. Full model vs Full model without finsqft
2. Full model vs Full model without bedroom
3. Full model vs Full model without lotisze

Each test is soundly rejected. The full model is much better with all three predictors.

drop1(m2, test = "F")

This test is also implemented in the Anova() function in the car package.

# install.packages("car")
library(car)
Anova(m2)

AIC and BIC

AIC (Akaike Information Criterion) is a statistic designed to help us choose a model with the best predictive power among a group of models. The value of AIC doesn’t really have any interpretation. It’s meant to be compared to other AIC values. The lower the AIC the better. We can use the AIC() function in R to compare multiple models. For example, mod_03 seems preferable to mod_02 because the AIC value is so much smaller.

AIC(mod_02, mod_03)

AIC is the log-likelihood of the model multiplied by -2 with 2 x df added to it. The 2 x df part is a penalty. “df” is short for “degrees of freedom” and is the number of parameters estimated in the model. This includes the residual standard error. For example, mod_03 has 5 degrees of freedom because there are 4 model coefficients and the residual standard error. We call this part a “penalty” because AIC can get bigger with more coefficients.

Log-likelihood is the log-transformed likelihood. Likelihood is the joint probability of the observed data as a function of the parameters of the chosen statistical model. Imagine turning the coefficients in the model like dials on a machine and trying to find the maximum likelihood. In other words, what combination of coefficients are most likely to produce the data we observed? R does not estimate linear model coefficients using maximum likelihood, but we can calculate the log likelihood after the fact using the logLik() function:

logLik(mod_03)

We can then calculate AIC “by-hand” as follows.

rbind("mod_02" = -2*logLik(mod_02) + 2*4, 
      "mod_03" = -2*logLik(mod_03) + 2*5)

The AIC is inclined to choose overly complex models, so some researchers prefer BIC (Bayesian Information Criterion), which places a bigger penalty on the number of predictors. Again use the BIC() function in the same way.

BIC(mod_02, mod_03)

The BIC is calculated the same as the AIC but with a different penalty, log(n), where n is the number of observations. Again we can calculate “by hand”:

rbind("mod_02" = -2*logLik(mod_02) + log(nrow(homes))*4, 
      "mod_03" = -2*logLik(mod_03) + log(nrow(homes))*5)

Of course, you don’t have to choose one. You can use both AIC and BIC and report both. They will often choose the same models.

Collinearity

When predictors are highly correlated (ie, have strong linear relationships), the precision of coefficient estimates can decline. This phenomenon is often referred to as collinearity or multicollinearity. Let’s demonstrate with a toy example. First we generate two variables, x1 and x2, that are highly correlated:

x1 <- seq(1,10,length = 100)
set.seed(123)
x2 <- 1 + 2*x1 + rnorm(100, sd = 0.2)
cor(x1, x2)  # calculate correlation; perfect correlation = 1

Now we generate a new variable, y, using x1 and x2 along with some noise and fit a linear model to recover the true coefficients of 2 and 3. Notice in the pairwise scatterplot that x1 and x2 both seem associated with y in the same way.

set.seed(321)
y <- 0.5 + 2*x1 + 3*x2 + rnorm(100, sd = 10)
d_collinear <- data.frame(y, x1, x2)
pairs(d_collinear)

When we fit the model, our estimated coefficients are way off and have enormous standard errors. The “true” values are 2 and 3. The model estimates are about 12 and -2! Recall we interpret the x1 coefficient as its effect on y holding x2 constant. But since x1 and x2 are highly correlated, it’s all but impossible to estimate the effect of x1 while holding x2 constant.

mod_colinear <- lm(y ~ x1 + x2)
summary(mod_colinear)

The most common way to check for and quantify collinearity after fitting a model is calculating variance inflation factors (VIF). The details of the calculation can be found with a web search, but the basic idea is that if one of the raw VIFs is greater than 10, then we may have evidence that collinearity is influencing coefficient estimates. Fortunately this is easy to do using the vif() function in the car package. The VIFs are sky high for our contrived predictors!

library(car)
vif(mod_colinear)

The square root of the VIF can be interpreted as how much larger the standard error of the coefficient is relative to similar uncorrelated data.

vif(mod_colinear) |> sqrt()

This says the standard error for both variables is about 29 times larger than it would have been without collinearity. The best solution in this case would be to simply drop one of the variables. It appears x1 is almost completely determined by x2 and vice versa. Knowing one means you know the other. In the extreme case, when two variables are perfectly correlated, the model fitting procedure will return NA for one of the variables, as demonstrated below. Notice the message: “(1 not defined because of singularities)”

x2 <- 2*x1  # perfect correlation
summary(lm(y ~ x1 + x2))

Let’s check our m2 model where we modeled log(totalprice) as a function of finsqft, bedroom and lotsize:

vif(m2)

These VIFs look very good. We might suspect that finsqft and number of bedrooms could be highly correlated, but the VIF checks out.

One approach to addressing collinearity concerns is to use a data reduction technique such as principal components analysis (PCA). When it works, this method basically takes several variables and reduces them to one or two summary scores. This may be preferable to arbitrarily dropping variables.

For models that are intended for making predictions, collinearity is not much of a concern.

Transformation guidelines

Above we log-transformed totalvalue to help meet modeling assumptions. Recall without the log transformation our residuals were large and skewed, which is a fancy way of saying our model was a bad fit. A good fitting model should have relatively small residuals with even scatter.

A log-transformation made sense for two reasons:

1. totalvalue was strictly positive, had a large upper bound, and covered several orders of magnitude.
2. changes in totalvalue according to the predictors were relative (multiplicative) and not absolute (additive), which corresponds to the natural log scale.

It’s important to note that not all skewed variables need to be transformed when it comes to linear modeling. The distributional assumptions are on the residuals. However there may be times when you need to investigate transformations other than the log when it comes to modeling. These almost always take the form of a power transformation (ie, raising your variable to a power using an exponent). Powers are usually symbolized with a Greek lambda (λ). A power of 0 translates to a log transformation.

Say your variable is y. A basic palette of possible power transformations include:

• λ = -1 1/y
• λ = -0.5 1/sqrt(y)
• λ = 0 log(y)
• λ = 0.5 sqrt(y)
• λ = 1 y (no transformation)
• λ = 2 y^2

The car function symbox() creates a visual assessment of which power makes the distribution reasonably symmetric. Below when we use it with totalvalue we see that the log transformation (λ = 0) does the best job of making the distribution more symmetric.

symbox(homes$totalvalue)

We can also use symbox() on a model object. For example, this produces essentially the same plot using the residuals of the model instead of totalvalue. Simply pipe the model into symbox().

lm(totalvalue ~ finsqft + bedroom + lotsize, data = homes) |>
  symbox()

A statistical “search” for the “best” power transformation can be performed with the powerTransform() function, also in the car package. The usual practice is to convert the result to the closest simple power listed above. For example, we can pipe the model result into powerTransform() and see the “best” transformation is about 0.16.

lm(totalvalue ~ finsqft + bedroom + lotsize, data = homes) |>
  powerTransform() 

0.16 is close to 0, so it makes sense to proceed with a log transformation. That greatly simplifies interpretation.

ggplot2 code for creating plot of simulated data

Here’s how to make the simulation plot using ggplot2. I find base R graphics easier for this type of plot.

sim2 <- simulate(m2, nsim = 50)
library(ggplot2)
library(tidyr)
sim2 %>% 
  pivot_longer(everything(), 
               names_to = "simulation", 
               values_to = "totalvalue") %>% 
  ggplot() +
  geom_density(mapping = aes(x = totalvalue, group = simulation),
               color = "grey80") + 
  geom_density(mapping = aes(x = log(totalvalue)),
               data = homes)  

How to make a function for simulating and plotting data for a linear model

base R

We basically copy and paste the original code in between the curly braces of the function() function. We call the function plot_sims but you can name it whatever you like. The changes are to the model name and number of simulations. We generalize those with arguments: mod and nsim. When we fit a model with lm, the data is stored with the model object by default and can be accessed as mod$model. The first column contains the model response, so we can access it with [,1] and use it to draw the density of the observed data.

plot_sims <- function(mod, nsim){
  sim <- simulate(mod, nsim = nsim)
  plot(density(mod$model[,1]))
  for(i in 1:nsim)lines(density(sim[[i]]), col = "grey80")
}
# try the function
plot_sims(mod = m2, nsim = 20)

ggplot2

Same idea as the previous function: we want to copy the original code in between the curly braces of the function() function. Except we now want to preface functions with their package name (eg, ggplot2::) so we can use the function without having the packages loaded. We also do away with the pipe %>% since it comes from yet another package (magrittr) but is accessible when tidyr is loaded. We extract the name of the response using resp <- names(mod$model)[1] and then use the use the .data pronoun (.data[[resp]]) from rlang package to use it with ggplot.

plot_sims <- function(mod, nsim){
  sim1 <- simulate(mod, nsim = nsim)
  resp <- names(mod$model)[1]
  sim1 <- tidyr::pivot_longer(sim1, everything(), 
               names_to = "simulation", 
               values_to = "totalvalue") 
  ggplot2::ggplot() +
  ggplot2::geom_density(mapping = ggplot2::aes(x = totalvalue, 
                                               group = simulation), 
                        color = "grey80", 
                        data = sim1) + 
  ggplot2::geom_density(mapping = ggplot2::aes(x = .data[[resp]]),
                        data = mod$model)  
}

plot_sims(m2, nsim = 65)

bayesplot package

The bayesplot package has a function for this called ppc_dens_overlay, but it’s a little weird to use for linear models because it was designed to be used with Bayesian models. However it’s not that hard to deploy. The first argument is simply the observed data. The second argument expects the simulations per row as opposed to per column. So we need to transpose, which we can do with the t() function. The result is a clean plot with the y-axis unlabeled since it really isn’t needed and a legend to distinguish between observed and simulated (or replicated) data.

# install.packages("bayesplot")
library(bayesplot)
ppc_dens_overlay(log(homes$totalvalue), t(sim2))

performance package

The easystats collection of packages “aims to provide a unifying and consistent framework to tame, discipline, and harness the scary R statistics and their pesky models.” https://easystats.github.io/easystats/ One of the packages is called {performance}, which provides the check_predictions() function for simulating data from a model and comparing it to the distribution of the original data. Two potential drawbacks at the time of this writing:

1. Can be slow for large data sets
2. log transformed data needs to be added as a variable to the data and modeled directly, as opposing to doing it on-the-fly as log(totalvalue)

Quick demo. Note the warning: “Minimum value of original data is not included in the replicated data.” This says the model is not generating data as small as the smallest value in the original data, which is 9600.

# install.packages("performance")
# install.packages("see")
library(performance)

# add log(totalvalue) to data and fit new model
homes$logtotalvalue <- log(homes$totalvalue)
m2a <- lm(logtotalvalue ~ finsqft + bedroom + lotsize, data = homes)
check_predictions(m2a)