Overview

Principle Component Analysis is widely used in data exploration, dimension reduction, data visualization. The aim is to transform original data into uncorrelated linear combinations of the original data while keeping the information contained in the data. High dimensional data tends to show clusters in lower dimensional view.

Clustering Analysis is another form of EDA. Here we are hoping to group data points which are close to each other within the groups and far away between different groups. Clustering using PC’s can be effective. Clustering analysis can be very subjective in the way we need to summarize the properties within each group.

Both PCA and Clustering Analysis are so called unsupervised learning. There is no response variables involved in the process.

For supervised learning, we try to find out how does a set of predictors relate to some response variable of the interest. Multiple regression is still by far, one of the most popular methods. We use a linear model as a working model for its simplicity and interpretability. It is important that we use domain knowledge as much as we can to determine the form of the response as well as the function format of the factors on the other hand.

Important Notice: This homework encompasses material from three modules. You will have a period of three weeks to complete it. Please manage your time accordingly.

0.1 Objectives

  • PCA
  • SVD
  • Clustering Analysis
  • Linear Regression

0.2 Review materials

  • Study Module 2: PCA
  • Study Module 3: Clustering Analysis
  • Study Module 4: Multiple regression (Including Simple regression as well)

0.3 Data needed

  • NLSY79.csv
  • brca_subtype.csv
  • brca_x_patient.csv

1 Case study 1: Self-esteem

Self-esteem generally describes a person’s overall sense of self-worthiness and personal value. It can play significant role in one’s motivation and success throughout the life. Factors that influence self-esteem can be inner thinking, health condition, age, life experiences etc. We will try to identify possible factors in our data that are related to the level of self-esteem.

In the well-cited National Longitudinal Study of Youth (NLSY79), it follows about 13,000 individuals and numerous individual-year information has been gathered through surveys. The survey data is open to public here. Among many variables we assembled a subset of variables including personal demographic variables in different years, household environment in 79, ASVAB test Scores in 81 and Self-Esteem scores in 81 and 87 respectively.

The data is store in NLSY79.csv.

Here are the description of variables:

Personal Demographic Variables

  • Gender: a factor with levels “female” and “male”
  • Education05: years of education completed by 2005
  • HeightFeet05, HeightInch05: height measurement. For example, a person of 5’10 will be recorded as HeightFeet05=5, HeightInch05=10.
  • Weight05: weight in lbs.
  • Income87, Income05: total annual income from wages and salary in 2005.
  • Job87 (missing), Job05: job type in 1987 and 2005, including Protective Service Occupations, Food Preparation and Serving Related Occupations, Cleaning and Building Service Occupations, Entertainment Attendants and Related Workers, Funeral Related Occupations, Personal Care and Service Workers, Sales and Related Workers, Office and Administrative Support Workers, Farming, Fishing and Forestry Occupations, Construction Trade and Extraction Workers, Installation, Maintenance and Repairs Workers, Production and Operating Workers, Food Preparation Occupations, Setters, Operators and Tenders, Transportation and Material Moving Workers

Household Environment

  • Imagazine: a variable taking on the value 1 if anyone in the respondent’s household regularly read magazines in 1979, otherwise 0
  • Inewspaper: a variable taking on the value 1 if anyone in the respondent’s household regularly read newspapers in 1979, otherwise 0
  • Ilibrary: a variable taking on the value 1 if anyone in the respondent’s household had a library card in 1979, otherwise 0
  • MotherEd: mother’s years of education
  • FatherEd: father’s years of education
  • FamilyIncome78

Variables Related to ASVAB test Scores in 1981

Test Description
AFQT percentile score on the AFQT intelligence test in 1981
Coding score on the Coding Speed test in 1981
Auto score on the Automotive and Shop test in 1981
Mechanic score on the Mechanic test in 1981
Elec score on the Electronics Information test in 1981
Science score on the General Science test in 1981
Math score on the Math test in 1981
Arith score on the Arithmetic Reasoning test in 1981
Word score on the Word Knowledge Test in 1981
Parag score on the Paragraph Comprehension test in 1981
Numer score on the Numerical Operations test in 1981

Self-Esteem test 81 and 87

We have two sets of self-esteem test, one in 1981 and the other in 1987. Each set has same 10 questions. They are labeled as Esteem81 and Esteem87 respectively followed by the question number. For example, Esteem81_1 is Esteem question 1 in 81.

The following 10 questions are answered as 1: strongly agree, 2: agree, 3: disagree, 4: strongly disagree

  • Esteem 1: “I am a person of worth”
  • Esteem 2: “I have a number of good qualities”
  • Esteem 3: “I am inclined to feel like a failure”
  • Esteem 4: “I do things as well as others”
  • Esteem 5: “I do not have much to be proud of”
  • Esteem 6: “I take a positive attitude towards myself and others”
  • Esteem 7: “I am satisfied with myself”
  • Esteem 8: “I wish I could have more respect for myself”
  • Esteem 9: “I feel useless at times”
  • Esteem 10: “I think I am no good at all”

1.1 Data preparation

Load the data. Do a quick EDA to get familiar with the data set. Pay attention to the unit of each variable. Are there any missing values?

This is a rich data set including information about 2,312 young adults over time (after cleaning). Originally, there are 46 columns with 2431 rows. Many columns are coded as integers with an exception for character values for Gender and Job05, and numerical for AFQT. We used a for-loop to looking for blank and missing values, and we found that there are blank values in column Job05. Besides that, there does not appear to be any missing values. There are integer value information of income, weight, height, years of education, scores on standardized tests, responses to a self confidence evaluation; as well as the categorical information of gender, job type, and access (yes or no) to magazines, newspapers, and libraries. The average respondent had 13.9 years of schooling, made $49,972 income, weighed 183 lbs, stood at 5.21 feet tall as of 2005, after scoring an average 55.3% on the AFQT exam and having parents with around 11 years of education each and family income of $21,547 in 1978. We see there were no obvious missing values so we proceed by using the unique function to find values that don’t make sense (negative income values). Additionally, we look at the range of values for the numerical categories to determine if values have the incorrect units, i.e. “63” on income is confusing when the other 99% of entries made at least $1000, making you think it was supposed to be 63k. We checked that all the yes/no questions only had “1” or “0” responses, that the AFQT percentage scores were between 0 and 100, that gender was only “female” or “male”, and that the only responses to the esteem questions were 1-4. Ultimately we remove entries with negative values or values that appear to be in the wrong units, adding up to only 4.9% of our original data set being lost to the cleaning process.

[1] “Empty values found in column Job05”

1.2 Self esteem evaluation

Let concentrate on Esteem scores evaluated in 87.

  1. First do a quick summary over all the Esteem variables. Pay attention to missing values, any peculiar numbers etc. How do you fix problems discovered if there is any? Briefly describe what you have done for the data preparation.

  2. Please note that higher scores on Esteem questions 1, 2, 4, 6, and 7 indicate higher self-esteem, whereas higher scores on the remaining questions suggest lower self-esteem. To maintain consistency, consider reversing the scores of certain Esteem questions. For example, if the esteem data is stored in data.esteem, you can use the code data.esteem[, c(1, 2, 4, 6, 7)] <- 5 - data.esteem[, c(1, 2, 4, 6, 7)] to invert the scores.

Data Preparation Summary: Among the data collected from the ‘87 Esteem survey, there were no missing values or peculiar numbers because each esteem question only had 1,2,3,4 as answers as expected. The only data preparation needed for the esteem variables was inverting scores on questions 1,2,4,6,7 because of the connotation being opposite of the other questions’ connotations. By doing this when we see an answer 4 to any of the questions that indicates lowest levels of self confidence.

  1. Write a brief summary with necessary plots about the 10 esteem measurements. Of all the responses from all subjects concerning all 10 esteem questions, .614% of responses were “Strongly Disagree”, 5.86% of responses were “Disagree”, 45.6% of responses were “Agree”, and 48% of responses were “Strongly Agree”. Additionally, the average score among all respondents and questions was 3.41, i.e. between the categories “Agree” and “Strongly Agree”. The average score for Questions 1,2,4,6,7 are less than 2 while Questions 3,5,8,9,10 are greater than 3. Questions associated with lower self-esteem have higher scores. This indicates that the majority of the time respondents indicate they either somewhat lack (“Agree”) or greatly lack (“Strongly Agree”) self confidence.

  1. Do esteem scores all positively correlated? Report the pairwise correlation table and write a brief summary.

All esteem scores are positively pairwise correlated. Overall question 9 has the lowest average correlation with all other questions and question 6 has the highest average correlation with all other questions. 1 and 2, 6 and 7, 4 and 2, 1 and 4, are >.5 while 3,5,9,10 are >.5 The highest correlation is between questions 1 and 2 (.7), and the lowest correlation is between questions 9 and 1 (.24).

  1. PCA on 10 esteem measurements. (centered but no scaling)
  1. Report the PC1 and PC2 loadings. Are they unit vectors? Are they orthogonal?
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Q1 0.232 -0.369 -0.065 0.021 0.394 0.049 -0.330 0.226 0.026 -0.698
Q2 0.242 -0.364 -0.059 -0.027 0.380 0.050 -0.301 0.238 -0.028 0.714
Q3 0.278 -0.150 -0.144 0.439 -0.108 -0.019 -0.285 -0.757 -0.140 0.020
Q4 0.256 -0.322 -0.048 -0.132 0.263 0.152 0.780 -0.196 -0.267 -0.027
Q5 0.310 -0.143 -0.093 0.567 -0.316 -0.436 0.266 0.424 0.111 0.011
Q6 0.313 -0.215 0.074 -0.358 -0.230 0.019 0.049 -0.196 0.792 0.015
Q7 0.301 -0.169 0.145 -0.502 -0.492 -0.242 -0.192 0.067 -0.514 -0.040
Q8 0.393 0.325 0.820 0.162 0.197 0.053 -0.002 -0.008 -0.017 -0.001
Q9 0.404 0.580 -0.429 -0.231 0.321 -0.396 0.002 -0.060 0.021 -0.006
Q10 0.378 0.261 -0.279 0.082 -0.289 0.750 -0.031 0.232 -0.058 -0.006

The PC1 and PC2 loadings are unit vectors and they are orthogonal.

  1. Are there good interpretations for PC1 and PC2? (If loadings are all negative, take the positive loadings for the ease of interpretation)

Loadings determine the contribution of each variable to the PCs. We see that for PC1 the loadings are all positive but are not all equivalent, but we can roughly say that that questions 1 and 2 (~.2 loading) almost equally contribute the least to PC1, questions 3-7 (~.3 loading) almost equally contribute the second least to PC1, and questions 8-10 (~.4 loading) almost equally contribute the most to PC1. As for PC2, we see that some loadings are positive and some are negative, meaning that PC2 is formed by a weighted average of the differences between scores on questions 1-7 and scores on questions 8-10. However, because the loadings are not equivalent in size we can not say PC2 is directly proportional to the difference in scores among these questions subsets. Generally, a higher PC1 score can be interpreted as higher scores on the esteem questions, while a higher PC2 score can be interpreted as higher scores on questions 8-9 and lower scores on questions 1-7.

  1. How is the PC1 score obtained for each subject? Write down the formula.

PC1 score for subject i = 0.232(subject i’s Q1 response)+ 0.242(subject i’s Q2 response) +0.278 (subject i’s Q3 response)+0.256(subject i’s Q4 response)+0.310(subject i’s Q5 response)+0.313(subject i’s Q6 response)+0.301(subject i’s Q7 response)+0.393(subject i’s Q8 response)+0.404(subject i’s Q9 response)+0.378(subject i’s Q10 response)

  1. Are PC1 scores and PC2 scores in the data uncorrelated? Yes! By definition they must be.
##       PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10
## PC1  1.67 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## PC2  0.00 0.46 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## PC3  0.00 0.00 0.32 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## PC4  0.00 0.00 0.00 0.27 0.00 0.00 0.00 0.00 0.00 0.00
## PC5  0.00 0.00 0.00 0.00 0.21 0.00 0.00 0.00 0.00 0.00
## PC6  0.00 0.00 0.00 0.00 0.00 0.19 0.00 0.00 0.00 0.00
## PC7  0.00 0.00 0.00 0.00 0.00 0.00 0.14 0.00 0.00 0.00
## PC8  0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.13 0.00 0.00
## PC9  0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.12 0.00
## PC10 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.07
  1. Plot PVE (Proportion of Variance Explained) and summarize the plot.

This plot shows the proportion of variance explained by each PC, which can help us decide how many PCs are needed. To make that decision, we look at the plot and apply the elbow rule, i.e. take the number of PCs when there is a sharp drop in the scree plot. Looking at this plot and applying the rule it seems we would only need the two leading PCs because PC3 is the elbow point.

  1. Also plot CPVE (Cumulative Proportion of Variance Explained). What proportion of the variance in the data is explained by the first two principal components?

About 60% of the variance in the data is explained by PC1 and PC2.

  1. PC’s provide us with a low dimensional view of the self-esteem scores. Use a biplot with the first two PC’s to display the data. Give an interpretation of PC1 and PC2 from the plot. (try ggbiplot if you could, much prettier!)

From the biplot, PC1 represents overall self-worth and positive attitude, while PC2 represents feelings of failure and dissatisfaction. The biplot provides a visual representation of how the original variables contribute to these principal components.

  1. Apply k-means to cluster subjects on the original esteem scores

  1. Find a reasonable number of clusters using within sum of squared with elbow rules. Two clusters appears to be a sufficient amount.

  2. Can you summarize common features within each cluster?

Those in Cluster 1 have lower average responses on all esteem questions, indicating they are more confident overall than those in Cluster 2 whom have higher average responses on the esteem survey.

  1. Can you visualize the clusters with somewhat clear boundaries? You may try different pairs of variables and different PC pairs of the esteem scores.

From testing all combinations of questions and comparing PC1 to PC2, the above relationships are where we can best see separation between the two clusters. Almost always we see blue with higher responses to questions. In the PC1 vs. PC2 plot we see blue with a higher PC1 score, which makes sense considering the higher responses on questions.

  1. We now try to find out what factors are related to self-esteem? PC1 of all the Esteem scores is a good variable to summarize one’s esteem scores. We take PC1 as our response variable.
  1. Prepare possible factors/variables:

  • EDA the data set first. Refer back to question 1.

    • Personal information: gender, education (05), log(income) in 87, job type in 87. One way to summarize one’s weight and height is via Body Mass Index which is defined as the body mass divided by the square of the body height, and is universally expressed in units of kg/m². Note, you need to create BMI first. Then may include it as one possible predictor.

  • Household environment: Imagazine, Inewspaper, Ilibrary, MotherEd, FatherEd, FamilyIncome78. Do set indicators Imagazine, Inewspaper and Ilibrary as factors.

    • You may use PC1 of ASVAB as level of intelligence
    1. Run a few regression models between PC1 of all the esteem scores and suitable variables listed in a). Find a final best model with your own clearly defined criterion.

  • How did you land this model? Run a model diagnosis to see if the linear model assumptions are reasonably met. We got to this model by starting with everything discussed as potential factors in self esteem above and parsing it down based on if those variables were reported to be significant. We also looked into different combinations of interactions but none of that proved fruitful.

## Warning: not plotting observations with leverage one:
##   668

We see that the homoskedasticity and linearity assumptions appear to be met since the points are somewhat randomly scattered in the residuals vs. fitted plot and the points follow the diagonal for the most part between the -2 to 2 quantiles in the QQ residuals plot. It is not perfect but this is much better than what we got from some of the other variable combinations we tried in models.

  • Write a summary of your findings. In particular, explain what and how the variables in the model affect one’s self-esteem.

In our model we find that PC1 of AFQT, years of education of the respondent as of 2005, jobs in the Health Care Technical and Support Occupations, Food Preparation and Serving Related Occupations, Cleaning and Building Service Occupations, and Transportation and Material Moving Workers are all significant to the 0 level, i.e. most significant. Additionally, being of the male gender; being in the fields of Engineers, Architects, Surveyers, Engineering and Related Technicians; Physical Scientists, Entertainment Attendants and Related Workers; Construction Trade and Extraction Workers, Installation, Maintenance and Repairs Workers; Setters, Operators and Tenders; and growing up in a household with access to a newspaper were all significant to the .01 level. Mother’s years of education; Farming, Fishing and Forestry Occupations; Sales and Related Workers; Teachers; and income in 1987 were all significant to the .05 level. To interpret the intercepts we need to keep in mind that a higher PC1 for the esteems questions indicates lower self confidence, so if a variable has a negative intercept that actually means it INCREASES self confidence. In our model results we see that being in the following fields INCREASE self confidence significantly: Engineers, Architects, Surveyers, Engineering and Related Technicians, Physical Scientists, Teachers, Health Care Technical and Support Occupations, Food Preparation and Serving Related Occupations, Cleaning and Building Service Occupations, and Transportation and Material Moving Workers, Entertainment Attendants and Related Workers; Construction Trade and Extraction Workers, Installation, Maintenance and Repairs Workers; Setters, Operators and Tenders.

2 Case Study 2: Breast cancer sub-type

The Cancer Genome Atlas (TCGA), a landmark cancer genomics program by National Cancer Institute (NCI), molecularly characterized over 20,000 primary cancer and matched normal samples spanning 33 cancer types. The genome data is open to public from the Genomic Data Commons Data Portal (GDC).

In this study, we focus on 4 sub-types of breast cancer (BRCA): basal-like (basal), Luminal A-like (lumA), Luminal B-like (lumB), HER2-enriched. The sub-type is based on PAM50, a clinical-grade luminal-basal classifier. (We had hoped to download the data for control groups for each type of the cancer. But failed to do so. Please let us know if you find the appropriate data.)

  • Luminal A cancers are low-grade, tend to grow slowly and have the best prognosis.
  • Luminal B cancers generally grow slightly faster than luminal A cancers and their prognosis is slightly worse.
  • HER2-enriched cancers tend to grow faster than luminal cancers and can have a worse prognosis, but they are often successfully treated with targeted therapies aimed at the HER2 protein.
  • Basal-like breast cancers or triple negative breast cancers do not have the three receptors that the other sub-types have so have fewer treatment options.

We will try to use mRNA expression data alone without the labels to classify 4 sub-types. Classification without labels or prediction without outcomes is called unsupervised learning. We will use K-means and spectrum clustering to cluster the mRNA data and see whether the sub-type can be separated through mRNA data.

We first read the data using data.table::fread() which is a faster way to read in big data than read.csv().

  1. Summary and transformation

  1. How many patients are there in each sub-type?

    ## brca_subtype
## Basal  Her2  LumA  LumB 
##   208    91   628   233

  2. Randomly pick 5 genes and plot the histogram by each sub-type.

    ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  3. Clean and transform the mRNA sequences by first remove gene with zero count and no variability and then apply logarithmic transform.

    # remove genes with 0 counts
sel_cols <- which(colSums(abs(brca)) != 0)
brca_sub <- brca[, sel_cols, with=F]
dim(brca_sub) # p was 19948
    ## [1]  1160 19669
    # use log scale to spread out all the cell numbers
brca_sub <- log2(as.matrix(brca_sub+1e-10))
  1. Apply kmeans on the transformed dataset with 4 centers (4 clusters) and output the discrepancy table between the real sub-type brca_subtype and the cluster labels.
1 2 3 4
Basal 17 3 1 187
Her2 9 23 43 16
LumA 86 147 395 0
LumB 22 105 104 2
  1. Spectrum clustering: to scale or not to scale?

  1. Apply PCA on the centered and scaled dataset. How many PCs should we use and why? You are encouraged to use irlba::irlba(). In order to do so please review the section about SVD in PCA module.

We should use 3 PCs because when we apply the elbow rule, there is a sharp drop in the scree plot after 3 PCs.

  1. Plot PC1 vs PC2 of the centered and scaled data and PC1 vs PC2 of the centered but unscaled data side by side. Should we scale or not scale for clustering process? Why? (Hint: to put plots side by side, use gridExtra::grid.arrange() or ggpubr::ggrrange() or egg::ggrrange() for ggplots; use fig.show="hold" as chunk option for base plots)

We should scale for clustering process because the kmeans clustering separates the breast cancer types in the scaled and centered data better than in the centered data, and it accounts for the distance.

  1. Spectrum clustering: center but do not scale the data

  1. Use the first 4 PCs of the centered and unscaled data and apply kmeans. Find a reasonable number of clusters using within sum of squared with the elbow rule.

Using within sum of squared with the elbow rule, it can be reasonable to use 4 clusters in this case.

  1. Choose an optimal cluster number and apply kmeans. Compare the real sub-type and the clustering label as follows: Plot scatter plot of PC1 vs PC2. Use point color to indicate the true cancer type and point shape to indicate the clustering label. Plot the kmeans centroids with black dots. Summarize how good is clustering results compared to the real sub-type.

Compared to the real sub-type, this clustering results appears to be mostly accurate with majority matching.

  1. Compare the clustering result from applying kmeans to the original data and the clustering result from applying kmeans to 4 PCs. Does PCA help in kmeans clustering? What might be the reasons if PCA helps?

PCA may help in kmeans clustering because it helps separate the clusters, but it may not be necessary.

  1. Now we have an x patient with breast cancer but with unknown sub-type. We have this patient’s mRNA sequencing data. Project this x patient to the space of PC1 and PC2. (Hint: remember we remove some gene with no counts or no variablity, take log and centered, then find its PC1 to PC4 scores) Plot this patient in the plot in b) with a black dot as well. Calculate the Euclidean distance between this patient and each of the centroid of the cluster. (Don’t forget the clusters are obtained by using 4 PC’s) Can you tell which sub-type this patient might have?

## [1] 288 695 260 508

This patient might have the Basal sub-type since it is located in/near that color cluster.

3 Case Study: Fuel Efficiency in Automobiles

Linda will refine this case study by the following Monday, Feb 12th)

What determines how fuel efficient a car is? Are Japanese cars more fuel efficient? To answer thes questions we will build various linear models using the Auto dataset from the book ISLR. The original dataset contains information for about 400 different cars built in various years. To get the data, first install the package ISLR which has been done in the first R-chunk. The Auto dataset should be loaded automatically. Original data source is here: https://archive.ics.uci.edu/ml/datasets/auto+mpg

Get familiar with this dataset first. Tip: you can use the command ?ISLR::Auto to view a description of the dataset. Our response variable will be MPG: miles per gallon.

3.1 EDA

This data set includes information about 392 vehicles, including mpg, cylinders, displacement, horsepower, weight, acceleration, model year, origin and car name. Among the cars in the sample, the miles per gallon ranges from 9 to 46.6 with a sample average of 23.4. The cars have between 3 to 8 cylinders, with an average of 5.47, referring to how powerful the car engines are with more cylinders indicating more power. The displacement ranges from 68 to 455 with mean 194, referring to the engine displacement in cubicle inches which is a measure of the engine’s size and indicates the total amount of air and fuel that the engine can draw in during one complete cycle. Horsepower ranges from 46 to 230 with an average of 104.5, and refers to a unit of power, i.e. the rate at which work is done. The weight of the cars range between 1613-5140 lbs with an average of 2978 lbs. The cars can accelerate from 0 to 60 mph in 8 to 24.8 seconds, with an average of 15.5 seconds. The cars’ model years range from 70 to 82. The origin of the car is a categorical variable we turn into a factor, where 1 indicates American, 2 indicates European, and 3 indicates Japanese. Lastly, there are vehicle names, which is also a categorical variable we turn into a factor.

In the pairwise plot, points that are red indicate American origin, blue indicate European origin, and green indicate Japanese origin. There are positive linear relationships between cylinders and displacement, cylinders and horsepower, cylinders and weight, displacement and horsepower, displacement and weight, and horsepower and weight. Miles per gallon decreases with increase in cylinders, displacement, horsepower, and weight. Time to acceleration also tends to decrease with higher displacement and horsepower. There are no real clear patterns surrouding year. American cars tend to have lower mpg and acceleration, but higher cylinders, displacement, horsepower, weight. European cars and Japanese cars are similar but opposite of American cars, with higher mpg and acceleration, lower cylinders, displacement, horsepower, weight.

  1. Explore the data, list the variables with clear definitions. Set each variable with its appropriate class. For example origin should be set as a factor. (see above for clear definitions)

  2. How many cars are included in this data set?

There are 392 cars in the data set.

  1. EDA, focus on pairwise plots and summary statistics. Briefly summarize your findings and any peculiarities in the data.

(see above for summary)

3.2 What effect does time have on MPG?

  1. Start with a simple regression of mpg vs. year and report R’s summary output. Is year a significant variable at the .05 level? State what effect year has on mpg, if any, according to this model.

Year is significant at the .05 level (actually up to 0 level). According to the model, on average, a one year increase in the model year of a car implies an average 1.23 increase in miles per gallon of a car.

  1. Add horsepower on top of the variable year to your linear model. Is year still a significant variable at the .05 level? Give a precise interpretation of the year’s effect found here. Year is still a significant variable at the .05 level (actually up to 0 level). Horsepower is also significant at the .05 level (actually up to 0 level). Here, we see that a one year increase in the model year of a car implies an average .65727 increase in miles per gallon of a car. We see that on average a one unit increase in horse power of a car indicates a decrease of -.13165 miles per gallon of a car.

  2. The two 95% CI’s for the coefficient of year differ among (a) and (b). How would you explain the difference to a non-statistician?

The 95% CI for year in model a is [1.06,1.4] and in model b is [.527,.788]. The interval for model a holds more values (is wider) than that of model b. This means in model a it is less clear what the true effect of the “year” variable on mpg is than in model b.

  1. Create a model with interaction by fitting lm(mpg ~ year * horsepower). Is the interaction effect significant at .05 level? Explain the year effect (if any).

The interaction effect is significant at the .05 level (actually up to 0 level), as well as year alone and horse power alone are significant. This indicates that the effect of year on mpg depends on the value of horsepower, or vice versa that the effect of horsepower on mpg depends on the year. We see that a one year increase in the model year of a car implies an average 2.19 increase in miles per gallon of a car when horse power is held constant.

3.3 Categorical predictors

Remember that the same variable can play different roles! Take a quick look at the variable cylinders, and try to use this variable in the following analyses wisely. We all agree that a larger number of cylinders will lower mpg. However, we can interpret cylinders as either a continuous (numeric) variable or a categorical variable.

  1. Fit a model that treats cylinders as a continuous/numeric variable. Is cylinders significant at the 0.01 level? What effect does cylinders play in this model?

Cylinders is significant at the .01 level (actually up to 0 level). Cylinders has a negative coefficient in this model, indicating that, on average, an increase in cylinders leads to a decrease in mpg for a car.

cylinders is already a continuous/numeric variable.

  1. Fit a model that treats cylinders as a categorical/factor. Is cylinders significant at the .01 level? What is the effect of cylinders in this model? Describe the cylinders effect over mpg.

3 cylinders (the intercept), 4 cylinders, and 8 cylinders are all significant at the .01 level. We see that cylinders 3, 4, and 5 have a positive effect on the mpg of the car on average, while cylinders 6 and 8 have a negative effect on the mpg of the car on average. It seems that being in one of the categories with less cylinders indicates the car has higher mpg, while categories with more cylinders indicate a car has lower mpg. Since categories with 5 and 6 cylinders are not significant and ‘on the edge’ this might be because they fall in the middle of two general car categories: fuel efficient and sports cars.

  1. What are the fundamental differences between treating cylinders as a continuous and categorical variable in your models? When cylinders is treated as a continuous variable, the model assumes that the effect of each additional cylinder is constant, so the estimated cylinders coefficient then represents the average change in the response variable for each one-unit increase in the number of cylinders. When cylinders is treated as categorical, the model estimates a coefficient for each amount of cylinders so these different cylinder categories are allowed to have different effects on the mpg of a car.

  2. Can you test the null hypothesis: fit0: mpg is linear in cylinders vs. fit1: mpg relates to cylinders as a categorical variable at .01 level?

The test is significant at the .01 level (down to the 0 level), indicating we can reject the null hypothesis and determine it is better for cylinders to be expressed as a categorical variable.

3.4 Results

Final modeling question: we want to explore the effects of each feature as best as possible. You may explore interactions, feature transformations, higher order terms, or other strategies within reason. The model(s) should be as parsimonious (simple) as possible unless the gain in accuracy is significant from your point of view.

  1. Describe the final model. Include diagnostic plots with particular focus on the model residuals and diagnoses.

The final model we developed uses year, weight, and horsepower as numerical variables and cylinders, origin, and brand as categorical variables. We develop the brands by taking just the first word of each car name. We achieve an R-squared value of .869, meaning 86.9% of the variance in mpg is explained by the model. Overall the model is significant.

  1. Summarize the effects found. The residuals vs. fitted plot passes the homoscedasticity check as the points are fairly equally varied across the horizontal axis, and this is further reassured by the randomness in the Scale-Location plot. In the residuals vs. fitted plot there might be slight concern for the linearity as the slight curve is showing up, but we judge this to not be major enough to throw out this model (plus all other models we tried were much less randomly dispersed). The Q-Q residuals plot gives us further reassurance to the goodness of our analysis given the data points follow the line within the -2 to 2 standard deviation range very well, suggesting it is normally distributed. The residuals vs. leverage plot is not perfect, but it does seem somewhat random and the suggested outliers are not repeated from other plots, so we accept it.

  2. Predict the mpg of the following car: A red car built in the US in 1983 that is 180 inches long, has eight cylinders, displaces 350 cu. inches, weighs 4000 pounds, and has a horsepower of 260. Also give a 95% CI for your prediction.

The red car with the given characteristics is predicted by our model to have 19 mpg, and the 95% CI of the prediction is from [11.7,26.2].

4 Simple Regression through simulations (Optional)

4.1 Linear model through simulations

This exercise is designed to help you understand the linear model using simulations. In this exercise, we will generate \((x_i, y_i)\) pairs so that all linear model assumptions are met.

Presume that \(\mathbf{x}\) and \(\mathbf{y}\) are linearly related with a normal error \(\boldsymbol{\varepsilon}\) , such that \(\mathbf{y} = 1 + 1.2\mathbf{x} + \boldsymbol{\varepsilon}\). The standard deviation of the error \(\varepsilon_i\) is \(\sigma = 2\).

We can create a sample input vector (\(n = 40\)) for \(\mathbf{x}\) with the following code:

# Generates a vector of size 40 with equally spaced values between 0 and 1, inclusive
x <- seq(0, 1, length = 40)

4.1.1 Generate data

Create a corresponding output vector for \(\mathbf{y}\) according to the equation given above. Use set.seed(1). Then, create a scatterplot with \((x_i, y_i)\) pairs. Base R plotting is acceptable, but if you can, please attempt to use ggplot2 to create the plot. Make sure to have clear labels and sensible titles on your plots.

4.1.2 Understand the model

  1. Find the LS estimates of \(\boldsymbol{\beta}_0\) and \(\boldsymbol{\beta}_1\), using the lm() function. What are the true values of \(\boldsymbol{\beta}_0\) and \(\boldsymbol{\beta}_1\)? Do the estimates look to be good?

Estimates aren’t bad, true values are 1 and 1.2

##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    1.331      0.557   2.390   0.0219
## x              0.906      0.959   0.946   0.3503
  1. What is your RSE for this linear model fit? Is it close to \(\sigma = 2\)?

1.79 is pretty close!

## [1] 1.79
  1. What is the 95% confidence interval for \(\boldsymbol{\beta}_1\)? Does this confidence interval capture the true \(\boldsymbol{\beta}_1\)?

[-1.03,2.85] Yes, it does capture 1.2

##  2.5 % 97.5 % 
##  -1.03   2.85
  1. Overlay the LS estimates and the true lines of the mean function onto a copy of the scatterplot you made above.

4.1.3 diagnoses

  1. Provide residual plot where fitted \(\mathbf{y}\)-values are on the x-axis and residuals are on the y-axis.

  1. Provide a normal QQ plot of the residuals.

  1. Comment on how well the model assumptions are met for the sample you used. Homoskedasticity and linearity checks appear to be met as the points are randomly disperes in the residuals vs. fitted plot and the points follow the diagonal line well in the QQ plot.

4.2 Understand sampling distribution and confidence intervals

This part aims to help you understand the notion of sampling statistics and confidence intervals. Let’s concentrate on estimating the slope only.

Generate 100 samples of size \(n = 40\), and estimate the slope coefficient from each sample.

  1. Summarize the LS estimates of \(\boldsymbol{\beta}_1\) (stored in results$b1). Does the sampling distribution agree with theory? Yes it does because 1.2 does end up in the 1st quartile to 3rd quartile range.

  2. How many of your 95% confidence intervals capture the true \(\boldsymbol{\beta}_1\)? Display your confidence intervals graphically.

96 of 100 of the confidence intervals captured the true \(\boldsymbol{\beta}_1\).

5 Statement on AI use

We used ChatGPT and Bard to help us figure out error messages and correct our code. Additionally, we used it to figure out how to put labels on bar charts and format within the Rmarkdown settings.