--- title: "Linear Regression" format: html number-sections: true theme: none html-math-method: mathml embed-resources: true from: markdown+lists_without_preceding_blankline --- ## Vector, list and dataframe We will be working with two data structures in R: **vector** and **list**. ### Vector We have seen how to create R vectors using the combine function `c`. ```{r} num_vec = c(1, 2, 3) ## a vector of numbers char_vec = c("I", "love", "econ") ## a vector of strings ``` In R, the vector is designed as a simple data structure: - All members of a vector must be of the same type (such as `char`, `numeric`, or `boolean`). - All vectors are flat. You cannot have a nested vector. ```{r} diff_types_vec = c(1, "a", TRUE) diff_types_vec ``` - R is very tolerant: instead of throwing an error, R silently converts both `1` (a numeric variable) and `TRUE` (a boolean variable) to characters. ```{r} diff_types_vec = c(1, TRUE) diff_types_vec ``` - In this case, R silently converts `TRUE` (a boolean variable) to `1` (a numeric variable). ```{r} nested_vec = c(1, c(10, 15), c(1,3, c(2,6))) nested_vec ``` - All vectors are flat. `c(1,c(3,4))` is the same as `c(1, 3, 4)`. ### List If you need an object containing variables of different types, use lists. ```{r} my_list = list(1, "a", TRUE) my_list ``` There are two ways to extract elements from a list: - `my_list[[1]]`: this returns a number `1` - `my_list[1]`: this returns a list `list(1)` Lists can be nested: ```{r} nested_list = list(1, "a", list(TRUE, 2)) nested_list ``` ### Lists vs vectors - When we are doing statistical computations, in most cases we only deal with numeric numbers. Use vectors. - Lists are more sophisticated than vectors. We use lists only when we want to store various kinds of information in a single object. - For example, in our example of linear regression last time, the object `my_lm` is indeed a list: - `my_lm = lm(y ~ x, data = my_data)` - Sometimes lists can be too powerful for our needs. It's common that we work with an Excel table of the following form: ``` id name age gender 1 Alice 20 F 2 Bob 23 M 3 Cindy 29 F ``` - Specifically, we have four vectors/variables in the above table: `id`, `name`, `age`, `gender`. R provides a special data structure for this kind of table-like object: dataframe. ### Dataframes `Dataframe` is a special kind of `list`. Like list, it can contain values of different kinds. We use `data.frame()` to create a dataframe: ```{r} df = data.frame( name = c("Alice", "Bob", "Cindy"), age = c(20, 23, 29), gender = c("F", "M", "F") ) df ``` To extract the `name` information in `df`, use `df$name` which works the same as `df[["name"]]`: ```{r} df$name ``` ## Boston dataset This example of linear regression is based on *Intro to Statistical Learning* (ISLR) ch3.6. It involves the **Boston** dataset, which pertains to Housing Values in Suburbs of Boston. The **Boston** data set is found in the `MASS` package. ```{r} library(MASS) ``` After loading `MASS`, the onject **Boston** is in your R workspace. Specifically, `Boston` is an object in R of type `list` (or `dataframe` here). ```{r} typeof(Boston) ``` ```{r} dim(Boston) ``` To get some feeling about the `Boston` dataset: - run `names(Boston)`, `head(Boston)`, `str(Boston)`, or `summary(Boston)` - run `?Boston` to see detailed documentation of the `Boston` dataset in `MASS` ```{r, eval=FALSE} ?Boston ``` ``` crim per capita crime rate by town. zn proportion of residential land zoned for lots over 25,000 sq.ft. indus proportion of non-retail business acres per town. chas Charles River dummy variable (= 1 if tract bounds river; 0 otherwise). nox nitrogen oxides concentration (parts per 10 million). rm average number of rooms per dwelling. age proportion of owner-occupied units built prior to 1940. dis weighted mean of distances to five Boston employment centres. rad index of accessibility to radial highways. tax full-value property-tax rate per $10,000. ptratio pupil-teacher ratio by town. black 1000*(Bk - 0.63)^2, where Bk is the proportion of blacks by town. lstat lower status of the population (percent). medv median value of owner-occupied homes in $1000s. ``` Another (much better) way to explore a dataset is to use plots. ### Exploring data using plots ```{r} plot(Boston$crim, Boston$medv, xlab="Crime rate", ylab="Median value of homes (in thousands)") ``` Based on the plot, it seems we should focus on those observations with `crim < 25` in plotting. ```{r,message=FALSE} library(dplyr) data_subset = tibble(crim=Boston$crim, medv=Boston$medv) data_subset = filter(data_subset, crim<25) ``` Here I use the `dplyr` package which provides some useful "verbs" (functions) for data manipulation. If you use the pipe operator `|>`, an equivalent (but more readable) way is: ```{r} df = as_tibble(Boston) |> select(crim, medv) |> filter(crim<25) ``` ```{r} plot(df$crim, df$medv, xlab="Crime rate (less than 25)", ylab="Median value of homes (in thousands)") ``` Plot each feature against `medv`: ```{r} library(reshape2) library(ggplot2) bosmelt <- melt(Boston, id="medv") ggplot(bosmelt, aes(x=value, y=medv))+ facet_wrap(~variable, scales="free")+ geom_point() ``` Lastly, we compute and visualize the correlation matrix: ```{r} corr_m = cor(Boston) round(corr_m, 2) ``` ```{r, message=FALSE} library(corrplot) ## visualizing the corr matrix corrplot(corr_m) ``` ### Linear regression We use the `lm()` function for running linear regression. For example, the syntax `lm(y ~ x1 + x2 + x3, data)` fits a model with three predictors: $y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3$. Here the `Boston` dataset contains 13 variables. You can manually specify all the 13 regressors. Fortunately, R provides the short-hand `.` for performing a regression using all of the predictors: ```{r} fit_all = lm(medv ~ ., data = Boston) ``` - As discussed above, `fit_all` is a list containing all the regression info. The ```{r} summary(fit_all) ``` Here the `summary` output still contains lots of info. We can use `summary(fit)$r.sq` to see the $R^2$. ```{r} summary(fit_all)$r.sq ## short for summary(fit_all)[["r.sq"]] ``` ### Excluding certain predictors From the summary, `age` and `indus` seem to have high p-values. We can run a regression excluding these two predictors: ```{r} fit_ex = lm(medv ~ . - age - indus, data=Boston) ``` ```{r} summary(fit_ex) ``` ### Interaction and nonlinear terms - The syntax `lstat:black` tells R to include the interaction term $\text{lstat}\times\text{black}$. - The syntax `lstat*age` simultaneously include `lstat`, `age`, and the interaction term `lstat:black` as predictors: ```{r} summary(lm(medv ~ lstat*age, data=Boston)) ``` To include the quadratic term `X^2`, use `I(X^2)`: - The special function `I()` is needed since the hat symbol `^` has a special meaning in a formula. ```{r} fit_nonlinear = lm(medv~lstat+I(lstat^2), data = Boston) summary(fit_nonlinear) ``` - The near-zero p-value associated with the quadratic term suggests that it leads to an improved model. ## Review questions 1. What's the difference between vector and list? Explain in your own words. 1. `m = matrix(1:4, nrow = 2)`. Is `m` a vector or a list? 1. Is the dataset `Boston` a vector or a list? Read about the dataset (`?Boston`) and answer the following questions. 4. How to know the number of rows and columns in dataset `Boston`? 1. Make some pairwise scatterplots of the predictors (columns) in this data set. Describe your findings. 1. Are any of the predictors associated with `crim` (per capita crime rate)? If so, explain the relationship. 1. Do any of the suburbs of Boston appear to have particularly high crime rates? Tax rates? Pupil-teacher ratios? Comment on the range of each predictor. 1. How many of the suburbs in this data set bound the Charles river?