Akshay Kumar Exam1 (Take Home)

Please hand this result in on Canvas no later than 11:59pm on Wednesday, March 12! Do not work in groups!

Consider the data from R package . We will use linear regression to investigate the relationship between variables in this data set and estimated performance (variable ). Do not use published performance as a predictor of performance in this problem.

# Load necessary libraries
library(ISLR2)

## Warning: package 'ISLR2' was built under R version 4.4.3

library(MASS)

## Warning: package 'MASS' was built under R version 4.4.3

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:ISLR2':
## 
##     Boston

library(tidyverse)

## Warning: package 'tidyr' was built under R version 4.4.2

## Warning: package 'readr' was built under R version 4.4.2

## Warning: package 'dplyr' was built under R version 4.4.2

## Warning: package 'lubridate' was built under R version 4.4.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(class)

## Warning: package 'class' was built under R version 4.4.3

library(car)

## Warning: package 'car' was built under R version 4.4.3

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.4.2

## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

library(e1071)

## Warning: package 'e1071' was built under R version 4.4.3

library(ggplot2)
library(dplyr)
library(TeachingDemos)

## Warning: package 'TeachingDemos' was built under R version 4.4.3

Linear Regression with cpus Data

Load and Prepare Data

data(cpus)

Investigate the relationship between variables in the dataset, both numerically and visually. Comment on the relationships you observe.

# Investigating relationships numerically and visually
summary(cpus)

##              name          syct             mmin            mmax      
##  ADVISOR 32/60 :  1   Min.   :  17.0   Min.   :   64   Min.   :   64  
##  AMDAHL 470/7A :  1   1st Qu.:  50.0   1st Qu.:  768   1st Qu.: 4000  
##  AMDAHL 470V/7 :  1   Median : 110.0   Median : 2000   Median : 8000  
##  AMDAHL 470V/7B:  1   Mean   : 203.8   Mean   : 2868   Mean   :11796  
##  AMDAHL 470V/7C:  1   3rd Qu.: 225.0   3rd Qu.: 4000   3rd Qu.:16000  
##  AMDAHL 470V/8 :  1   Max.   :1500.0   Max.   :32000   Max.   :64000  
##  (Other)       :203                                                   
##       cach            chmin            chmax             perf       
##  Min.   :  0.00   Min.   : 0.000   Min.   :  0.00   Min.   :   6.0  
##  1st Qu.:  0.00   1st Qu.: 1.000   1st Qu.:  5.00   1st Qu.:  27.0  
##  Median :  8.00   Median : 2.000   Median :  8.00   Median :  50.0  
##  Mean   : 25.21   Mean   : 4.699   Mean   : 18.27   Mean   : 105.6  
##  3rd Qu.: 32.00   3rd Qu.: 6.000   3rd Qu.: 24.00   3rd Qu.: 113.0  
##  Max.   :256.00   Max.   :52.000   Max.   :176.00   Max.   :1150.0  
##                                                                     
##     estperf       
##  Min.   :  15.00  
##  1st Qu.:  28.00  
##  Median :  45.00  
##  Mean   :  99.33  
##  3rd Qu.: 101.00  
##  Max.   :1238.00  
##

pairs(cpus[, -c(1, 8)]) # excluding 'name' and 'perf'

# Correlation matrix
cor(cpus[, -c(1, 8)])

##               syct       mmin       mmax       cach      chmin      chmax
## syct     1.0000000 -0.3356422 -0.3785606 -0.3209998 -0.3010897 -0.2505023
## mmin    -0.3356422  1.0000000  0.7581573  0.5347291  0.5171892  0.2669074
## mmax    -0.3785606  0.7581573  1.0000000  0.5379898  0.5605134  0.5272462
## cach    -0.3209998  0.5347291  0.5379898  1.0000000  0.5822455  0.4878458
## chmin   -0.3010897  0.5171892  0.5605134  0.5822455  1.0000000  0.5482812
## chmax   -0.2505023  0.2669074  0.5272462  0.4878458  0.5482812  1.0000000
## estperf -0.2883956  0.8192915  0.9012024  0.6486203  0.6105802  0.5921556
##            estperf
## syct    -0.2883956
## mmin     0.8192915
## mmax     0.9012024
## cach     0.6486203
## chmin    0.6105802
## chmax    0.5921556
## estperf  1.0000000

# Explanation:
# The summary provides basic statistics. The pairs plot visually shows relationships among all predictors and estperf.
# The correlation matrix numerically quantifies linear relationships.

Use either methods commonly used in the book/lecture notes to build a linear regression model predicting estimated performance from predictors in the dataset. Do not consider in this modeling approach. Explain the process used to arrive at your final model.

#  Building a linear regression model predicting estperf
initial_model <- lm(estperf ~ syct + mmin + mmax + cach + chmin + chmax, data = cpus)
stepwise_model <- step(initial_model <- lm(estperf ~ ., data = cpus[, -c(1,8)]), direction = "both")

## Start:  AIC=1615.2
## estperf ~ syct + mmin + mmax + cach + chmin + chmax
## 
##         Df Sum of Sq    RSS    AIC
## - chmin  1       146 444155 1613.3
## <none>               444009 1615.2
## - cach   1     45173 489182 1633.5
## - syct   1     51014 495022 1635.9
## - chmax  1    107238 551246 1658.4
## - mmin   1    220724 664732 1697.5
## - mmax   1    379414 823423 1742.3
## 
## Step:  AIC=1613.27
## estperf ~ syct + mmin + mmax + cach + chmax
## 
##         Df Sum of Sq    RSS    AIC
## <none>               444155 1613.3
## + chmin  1       146 444009 1615.2
## - cach   1     47080 491235 1632.3
## - syct   1     51380 495535 1634.2
## - chmax  1    117044 561199 1660.2
## - mmin   1    227105 671260 1697.6
## - mmax   1    379580 823735 1740.4

summary(stepwise_model)

## 
## Call:
## lm(formula = estperf ~ syct + mmin + mmax + cach + chmax, data = cpus[, 
##     -c(1, 8)])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -136.633  -23.149    1.574   22.957  287.869 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.660e+01  6.257e+00 -10.643  < 2e-16 ***
## syct         6.613e-02  1.365e-02   4.846 2.50e-06 ***
## mmin         1.424e-02  1.397e-03  10.188  < 2e-16 ***
## mmax         6.584e-03  4.999e-04  13.171  < 2e-16 ***
## cach         4.871e-01  1.050e-01   4.639 6.28e-06 ***
## chmax        1.187e+00  1.623e-01   7.314 5.88e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 46.78 on 203 degrees of freedom
## Multiple R-squared:  0.9108, Adjusted R-squared:  0.9086 
## F-statistic: 414.8 on 5 and 203 DF,  p-value: < 2.2e-16

# Explanation:
# We build an initial model with all predictors (excluding name and perf), then use stepwise regression to select the best subset of predictors based on AIC.

\vspace{20pt}

Create a residual plot using this model and comment on it’s features. Do any of the assumptions of linear regression seem to be violated? What might be done to adjust our model? Adjust the model if necessary by considering various residual plots, updating the model, and assessing residual plots using the updated model.

# Residual plot and checking assumptions
plot(stepwise_model, which = 1) # residual vs fitted values plot

# Explanation:
# The residual plot helps assess assumptions like homoscedasticity and linearity. If patterns or heteroscedasticity are present, transformations may be needed.

# Adjusting model if necessary (log-transforming response if residuals show issues)
log_model <- lm(log(estperf) ~ ., data = cpus[, -c(1,8)])
plot(log_model$fitted.values, log_model$residuals, main="Residual Plot", xlab="Fitted", ylab="Residuals")

# Explanation:
# Log transformation can stabilize variance if residuals indicate heteroscedasticity or nonlinearity.

\vspace{20pt}

How well does the final model fit the data? Comment on some model fit criteria from the model built in c)

# Model fit criteria
summary(log_model)$r.squared # R-squared

## [1] 0.9470454

AIC(log_model)

## [1] -33.25843

# Explanation:
# R-squared and adjusted R-squared values indicate how well the model fits the data. Higher values indicate better fit.

\vspace{20pt}

Interpret all variables in your final model using complete sentences, making sure to account for the fact that this may be a multivariable model. Give interpretations in terms of as meaningful of units as possible (it may not be possible to use seconds for cycle time - the answer is too large, but you may use MB instead of kB, for instance). Adjust interpretations as needed, both for units, and the fact that our outcome has been log transformed (how do we get to the raw data values from a log transformation? Start by thinking: what is the inverse of the log function???)

# Build the final model (assuming stepwise regression was used earlier)
final_model <- lm(log(estperf) ~ syct + mmin + mmax + cach + chmin + chmax, data = cpus)

# Extract coefficients
coefficients <- coef(final_model)

# Exponentiate coefficients to interpret them on the original scale of estperf
exp_coefficients <- exp(coefficients)

# Display coefficients and their exponentiated values
data.frame(
  Variable = names(coefficients),
  Coefficient = coefficients,
  Exp_Coefficient = exp_coefficients
)

Model Coefficients:
- We use coef() to extract coefficients from the final model. These coefficients are on the log scale.
Exponentiation:
- Using exp(), we transform the coefficients back to their original scale. This allows us to interpret them as multiplicative effects on estperf.
Interpretation:
- For each variable (e.g., syct, mmin, etc.), we calculate how a one-unit increase affects estperf multiplicatively.
- For variables like mmin and cach, we convert units (e.g., KB to MB or increments of 100 KB) for more meaningful interpretations.
  
  \vspace{20pt}

Calculate indices that help assess multicollinearity between predictors in your final model. Is there evidence of multicollinearity? What does this imply, and should you take action? Take action if appropriate.

# Fit the final regression model (log-transformed response)
final_model <- lm(log(estperf) ~ syct + mmin + mmax + cach + chmin + chmax, data = cpus)

# Calculate Variance Inflation Factor (VIF)
vif_values <- vif(final_model)
print(vif_values)

##     syct     mmin     mmax     cach    chmin    chmax 
## 1.201644 2.902002 3.274002 1.858349 1.966234 1.891392

# Explanation:
# The `vif()` function calculates Variance Inflation Factors for each predictor in the model.
# VIF values greater than 5 or 10 indicate potential multicollinearity issues.

\vspace{20pt}

Are there any outliers or influential observations in this data? Calculate relevant indices or provide visualizations to justify your answer. Make sure to use rules of thumb discussed in class if necessary for interpretations.

# Fit the final regression model (log-transformed response)
final_model <- lm(log(estperf) ~ syct + mmin + mmax + cach + chmin + chmax, data = cpus)

# Calculate Cook's distance
cooks_d <- cooks.distance(final_model)

# Plot Cook's distance
plot(cooks_d, type="h", main="Cook's Distance", ylab="Cook's Distance")
abline(h = 4/(nrow(cpus)-length(final_model$coefficients)-1), col="red", lty=2)
text(x=1:length(cooks_d), y=cooks_d, labels=ifelse(cooks_d > 4/(nrow(cpus)-length(final_model$coefficients)-1), names(cooks_d), ""), pos=4, cex=0.7)

# Identify influential points based on Cook's distance
influential_points <- which(cooks_d > (4/(nrow(cpus)-length(final_model$coefficients)-1)))
print(influential_points)

##   9  10  20  83  99 123 124 138 157 199 200 
##   9  10  20  83  99 123 124 138 157 199 200

# Calculate leverage values
leverage_values <- hatvalues(final_model)
leverage_cutoff <- 3*(length(final_model$coefficients))/nrow(cpus)

# Plot leverage values
plot(leverage_values, type="h", main="Leverage Values", ylab="Leverage")
abline(h = leverage_cutoff, col="blue", lty=2)
text(x=1:length(leverage_values), y=leverage_values, labels=ifelse(leverage_values > leverage_cutoff, names(leverage_values), ""), pos=4, cex=0.7)

# Identify high-leverage points
high_leverage_points <- which(leverage_values > leverage_cutoff)
print(high_leverage_points)

##   1   9  10  83 123 124 138 157 196 197 198 199 200 
##   1   9  10  83 123 124 138 157 196 197 198 199 200

# Studentized residuals to detect outliers
student_resid <- rstudent(final_model)

# Plot studentized residuals
plot(student_resid, type="h", main="Studentized Residuals", ylab="Studentized Residual")
abline(h=c(-3,3), col="purple", lty=2)
text(x=1:length(student_resid), y=student_resid, labels=ifelse(abs(student_resid)>3,names(student_resid),""), pos=4, cex=0.7)

# Identify outliers based on studentized residuals (> |3|)
outliers <- which(abs(student_resid) > 3)
print(outliers)

##   9  10 157 200 
##   9  10 157 200

Cook’s Distance: Measures influence of observations on regression coefficients. We use a common cutoff value 4/(n−k−1)4/(n−k−1) (where nn is number of observations and kk is number of predictors). Observations above this line are influential.
Leverage Values: Measure how far an observation’s predictor values are from the mean predictor values. Observations with leverage greater than 3(k+1)/n3(k+1)/n are considered high-leverage points.
Studentized Residuals: Identify outliers in terms of vertical distance from the fitted regression line. Observations with absolute studentized residuals greater than 3 are considered outliers.

\vspace{40pt}

Consider the data from R package . We will investigate the relationship between low birthweight and the predictors in the data using logistic regression and discriminant analysis.

Investigate the relationship between variables in the dataset. Do you see anything surprising? Use both numeric and visual summaries. Create and comment on visualizations specifically between the outcome variable and predictor/independent variables. Also, notice that qualitative/categorical variables should be visualized in an alternative manner, not just scatterplots/correlations as in the case of quantitative variables.

# Load the birthwt dataset
data(birthwt)

# a) Investigate relationships numerically and visually

# Numeric summaries
summary(birthwt)

##       low              age             lwt             race      
##  Min.   :0.0000   Min.   :14.00   Min.   : 80.0   Min.   :1.000  
##  1st Qu.:0.0000   1st Qu.:19.00   1st Qu.:110.0   1st Qu.:1.000  
##  Median :0.0000   Median :23.00   Median :121.0   Median :1.000  
##  Mean   :0.3122   Mean   :23.24   Mean   :129.8   Mean   :1.847  
##  3rd Qu.:1.0000   3rd Qu.:26.00   3rd Qu.:140.0   3rd Qu.:3.000  
##  Max.   :1.0000   Max.   :45.00   Max.   :250.0   Max.   :3.000  
##      smoke             ptl               ht                ui        
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.0000  
##  Median :0.0000   Median :0.0000   Median :0.00000   Median :0.0000  
##  Mean   :0.3915   Mean   :0.1958   Mean   :0.06349   Mean   :0.1481  
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:0.0000  
##  Max.   :1.0000   Max.   :3.0000   Max.   :1.00000   Max.   :1.0000  
##       ftv              bwt      
##  Min.   :0.0000   Min.   : 709  
##  1st Qu.:0.0000   1st Qu.:2414  
##  Median :0.0000   Median :2977  
##  Mean   :0.7937   Mean   :2945  
##  3rd Qu.:1.0000   3rd Qu.:3487  
##  Max.   :6.0000   Max.   :4990

# Correlation matrix for quantitative variables
quant_vars <- birthwt[, c("age", "lwt", "bwt")] # Select quantitative variables
cor_matrix <- cor(quant_vars)
print(cor_matrix)

##            age       lwt        bwt
## age 1.00000000 0.1800732 0.09031781
## lwt 0.18007315 1.0000000 0.18573328
## bwt 0.09031781 0.1857333 1.00000000

# Visualizations for quantitative variables
pairs(quant_vars, main = "Scatterplot Matrix of Quantitative Variables")

# Boxplot for outcome variable (low) vs age
ggplot(birthwt, aes(x = factor(low), y = age)) +
  geom_boxplot() +
  labs(x = "Low Birthweight (0 = No, 1 = Yes)", y = "Mother's Age") +
  ggtitle("Boxplot of Age by Low Birthweight")

# Boxplot for outcome variable (low) vs mother's weight at last menstrual period (lwt)
ggplot(birthwt, aes(x = factor(low), y = lwt)) +
  geom_boxplot() +
  labs(x = "Low Birthweight (0 = No, 1 = Yes)", y = "Mother's Weight at Last Menstrual Period (lwt)") +
  ggtitle("Boxplot of Weight by Low Birthweight")

# Barplots for categorical predictors
categorical_vars <- c("race", "smoke", "ht", "ui")
for (var in categorical_vars) {
  print(
    ggplot(birthwt, aes_string(x = var, fill = "factor(low)")) +
      geom_bar(position = "dodge") +
      labs(x = var, fill = "Low Birthweight") +
      ggtitle(paste("Barplot of", var, "by Low Birthweight")) +
      theme_minimal()
  )
}

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Explanation:
# - `summary()` provides numeric summaries for all variables.
# - `cor()` computes correlations between quantitative variables.
# - `pairs()` creates scatterplots to explore relationships between quantitative variables.
# - Boxplots show how continuous predictors (e.g., age, lwt) vary with the outcome (`low`).
# - Barplots visualize relationships between categorical predictors (e.g., race, smoke) and `low`.

Numeric Summaries:
- Use summary() to summarize all variables.
- Compute correlations for quantitative variables to understand linear relationships.
Visual Summaries:
- Scatterplots (pairs) are used for quantitative variables.
- Boxplots are used to compare continuous predictors across levels of the binary outcome (low).
- Barplots are used to show distributions of categorical predictors like race, smoke, etc., across levels of low.
  
  \vspace{20pt}

Fit a logistic regression model using methods discussed in class/the book, similar to as in problem 1). Be careful to understand each variable in to avoid including variables that are not logically acceptable for inclusion in the model.

# Convert categorical variables to factors for logistic regression
birthwt$race <- factor(birthwt$race, labels = c("White", "Black", "Other"))
birthwt$smoke <- factor(birthwt$smoke, labels = c("Non-Smoker", "Smoker"))
birthwt$ht <- factor(birthwt$ht, labels = c("No", "Yes"))
birthwt$ui <- factor(birthwt$ui, labels = c("No", "Yes"))
birthwt$low <- factor(birthwt$low, labels = c("No", "Yes"))

# Fit an initial logistic regression model with all predictors
logit_full <- glm(low ~ age + lwt + race + smoke + ptl + ht + ui + ftv, 
                  data = birthwt, family = binomial())

# Summary of the full model
summary(logit_full)

## 
## Call:
## glm(formula = low ~ age + lwt + race + smoke + ptl + ht + ui + 
##     ftv, family = binomial(), data = birthwt)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.480623   1.196888   0.402  0.68801   
## age         -0.029549   0.037031  -0.798  0.42489   
## lwt         -0.015424   0.006919  -2.229  0.02580 * 
## raceBlack    1.272260   0.527357   2.413  0.01584 * 
## raceOther    0.880496   0.440778   1.998  0.04576 * 
## smokeSmoker  0.938846   0.402147   2.335  0.01957 * 
## ptl          0.543337   0.345403   1.573  0.11571   
## htYes        1.863303   0.697533   2.671  0.00756 **
## uiYes        0.767648   0.459318   1.671  0.09467 . 
## ftv          0.065302   0.172394   0.379  0.70484   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 201.28  on 179  degrees of freedom
## AIC: 221.28
## 
## Number of Fisher Scoring iterations: 4

# Explanation:
# - The response variable `low` (low birthweight) is binary and modeled using logistic regression.
# - Predictors include both continuous variables (e.g., age, lwt) and categorical variables (e.g., race, smoke).
# - Categorical variables are converted to factors to ensure proper treatment in the model.

# Perform stepwise selection to refine the model
logit_final <- step(logit_full, direction = "both")

## Start:  AIC=221.28
## low ~ age + lwt + race + smoke + ptl + ht + ui + ftv
## 
##         Df Deviance    AIC
## - ftv    1   201.43 219.43
## - age    1   201.93 219.93
## <none>       201.28 221.28
## - ptl    1   203.83 221.83
## - ui     1   204.03 222.03
## - race   2   208.75 224.75
## - lwt    1   206.80 224.80
## - smoke  1   206.91 224.91
## - ht     1   208.81 226.81
## 
## Step:  AIC=219.43
## low ~ age + lwt + race + smoke + ptl + ht + ui
## 
##         Df Deviance    AIC
## - age    1   201.99 217.99
## <none>       201.43 219.43
## - ptl    1   203.95 219.95
## - ui     1   204.11 220.11
## + ftv    1   201.28 221.28
## - race   2   208.77 222.77
## - lwt    1   206.81 222.81
## - smoke  1   206.92 222.92
## - ht     1   208.81 224.81
## 
## Step:  AIC=217.99
## low ~ lwt + race + smoke + ptl + ht + ui
## 
##         Df Deviance    AIC
## <none>       201.99 217.99
## - ptl    1   204.22 218.22
## - ui     1   204.90 218.90
## + age    1   201.43 219.43
## + ftv    1   201.93 219.93
## - smoke  1   207.73 221.73
## - lwt    1   208.11 222.11
## - race   2   210.31 222.31
## - ht     1   209.46 223.46

# Summary of the final model after stepwise selection
summary(logit_final)

## 
## Call:
## glm(formula = low ~ lwt + race + smoke + ptl + ht + ui, family = binomial(), 
##     data = birthwt)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.086550   0.951760  -0.091  0.92754   
## lwt         -0.015905   0.006855  -2.320  0.02033 * 
## raceBlack    1.325719   0.522243   2.539  0.01113 * 
## raceOther    0.897078   0.433881   2.068  0.03868 * 
## smokeSmoker  0.938727   0.398717   2.354  0.01855 * 
## ptl          0.503215   0.341231   1.475  0.14029   
## htYes        1.855042   0.695118   2.669  0.00762 **
## uiYes        0.785698   0.456441   1.721  0.08519 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 201.99  on 181  degrees of freedom
## AIC: 217.99
## 
## Number of Fisher Scoring iterations: 4

# Explanation:
# - The `step()` function performs stepwise model selection based on AIC to identify the most relevant predictors.
# - The final model includes only those predictors that contribute significantly to explaining the variability in low birthweight.

# Check coefficients and interpret odds ratios
exp(coef(logit_final))

## (Intercept)         lwt   raceBlack   raceOther smokeSmoker         ptl 
##   0.9170901   0.9842205   3.7648926   2.4524265   2.5567241   1.6540303 
##       htYes       uiYes 
##   6.3919640   2.1939368

# Explanation:
# - The `exp()` function is used to calculate odds ratios from the logistic regression coefficients.
# - Odds ratios represent the multiplicative change in odds of low birthweight for a one-unit increase in a predictor (or change in category for categorical predictors).

Logistic Regression:
- The binary outcome variable low is modeled using logistic regression (glm() with family = bionomial())
- Predictors include both continuous variables (e.g., age, lwt) and categorical variables (e.g., race, smoke, ht, ui).
Stepwise Selection:
- Stepwise selection (step()) refines the model by iteratively adding or removing predictors based on AIC.
- This ensures that only statistically significant and relevant predictors are retained in the final model.
Odds Ratios:
- Exponentiating coefficients (exp(coef())) provides odds ratios, which are easier to interpret than raw coefficients.

Output:

Initial Model: Includes all predictors.
Final Model: Reduced set of predictors after stepwise selection.
Odds Ratios: Show how each predictor affects the odds of low birthweight.
\vspace{20pt}

What do you notice regarding the variables and . What is your logistic regression model in b) (perhaps before performing variable selection) implicitly assuming regarding these variables’ effects on the log odds of giving birth to a low weight baby? Are these assumptions realistic?

# Convert categorical variables to factors for logistic regression
birthwt$race <- factor(birthwt$race, labels = c("White", "Black", "Other"))
birthwt$smoke <- factor(birthwt$smoke, labels = c("Non-Smoker", "Smoker"))
birthwt$ht <- factor(birthwt$ht, labels = c("No", "Yes"))
birthwt$ui <- factor(birthwt$ui, labels = c("No", "Yes"))
birthwt$low <- factor(birthwt$low, labels = c("No", "Yes"))

# Fit a logistic regression model including ptl and ftv
logit_model <- glm(low ~ age + lwt + race + smoke + ptl + ht + ui + ftv, 
                   data = birthwt, family = binomial())

# Summary of the model
summary(logit_model)

## 
## Call:
## glm(formula = low ~ age + lwt + race + smoke + ptl + ht + ui + 
##     ftv, family = binomial(), data = birthwt)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.480623   1.196888   0.402  0.68801   
## age         -0.029549   0.037031  -0.798  0.42489   
## lwt         -0.015424   0.006919  -2.229  0.02580 * 
## raceBlack    1.272260   0.527357   2.413  0.01584 * 
## raceOther    0.880496   0.440778   1.998  0.04576 * 
## smokeSmoker  0.938846   0.402147   2.335  0.01957 * 
## ptl          0.543337   0.345403   1.573  0.11571   
## htYes        1.863303   0.697533   2.671  0.00756 **
## uiYes        0.767648   0.459318   1.671  0.09467 . 
## ftv          0.065302   0.172394   0.379  0.70484   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 201.28  on 179  degrees of freedom
## AIC: 221.28
## 
## Number of Fisher Scoring iterations: 4

# Check distribution of ptl and ftv
table(birthwt$ptl)

## 
##   0   1   2   3 
## 159  24   5   1

table(birthwt$ftv)

## 
##   0   1   2   3   4   6 
## 100  47  30   7   4   1

# Explanation:
# - The `summary()` function provides coefficients for all predictors, including `ptl` and `ftv`.
# - The `table()` function shows the distribution of values for `ptl` and `ftv` in the dataset.

# Visualize relationships between ptl/ftv and low birthweight
library(ggplot2)

ggplot(birthwt, aes(x = factor(ptl), fill = low)) +
  geom_bar(position = "dodge") +
  labs(x = "Number of Previous Preterm Labors (ptl)", fill = "Low Birthweight") +
  ggtitle("Barplot of ptl by Low Birthweight")

ggplot(birthwt, aes(x = factor(ftv), fill = low)) +
  geom_bar(position = "dodge") +
  labs(x = "Number of Physician Visits in First Trimester (ftv)", fill = "Low Birthweight") +
  ggtitle("Barplot of ftv by Low Birthweight")

# Explanation:
# - Barplots show how the distribution of `ptl` and `ftv` varies by low birthweight status.

Key Observations:

Implicit Assumptions:
- In logistic regression, continuous or ordinal variables like ptl and ftv are assumed to have a linear effect on the log odds of the outcome (low).
- For example, each unit increase in ptl or ftv is assumed to have a constant multiplicative effect on the odds of low birthweight.
Realism of Assumptions:
- The distribution of ptl and ftv may not support this assumption. For instance:
- If most values for ptl are zero or one, treating it as a continuous variable may not capture its true effect.
- Similarly, if ftv has many zeros or small values, its effect may not be linear.
  
  \vspace{20pt}

Create a new variable for named which is more useful for analysis. Keep in mind that with very small sample sizes, it may be worthwhile to collapse multiple categories.

# Create a new variable ptl2
# Collapse categories: if ptl > 0, set to '1+', otherwise '0'
birthwt$ptl2 <- factor(ifelse(birthwt$ptl > 0, "1+", "0"))

# Check the distribution of the new variable ptl2
table(birthwt$ptl2)

## 
##   0  1+ 
## 159  30

# Explanation:
# - The `ifelse()` function is used to create a binary variable `ptl2`.
#   - If `ptl > 0`, it is categorized as "1+" (indicating one or more previous preterm labors).
#   - Otherwise, it is categorized as "0" (indicating no previous preterm labors).
# - The `factor()` function ensures that the new variable is treated as a categorical variable.
# - The `table()` function provides the frequency distribution of the new variable to verify its creation.

Key Points:

Purpose:
- Collapsing categories reduces the number of levels in ptl, which is useful when dealing with small sample sizes.
- This binary categorization simplifies analysis while retaining meaningful distinctions.
Output:
- The table() function will display the counts for each category (“0” and “1+”) in ptl2.
  
  \vspace{20pt}
e. Create a new variable for named which is more useful for analysis. Keep in mind that with very small sample sizes, it may be worthwhile to collapse multiple categories. Also, it may be helpful to form tables which summarize low birthweight probabilities by levels of the variable in order to better understand the relationship between probability of low birthweight and the newly created variable.

# Map 'low' levels to numeric values: "No" -> 0, "Yes" -> 1
birthwt$low_numeric <- ifelse(birthwt$low == "No", 0, ifelse(birthwt$low == "Yes", 1, NA))

# Check for any NAs introduced during conversion
if (any(is.na(birthwt$low_numeric))) {
  cat("Warning: NAs were introduced during conversion of 'low' to numeric.\n")
  print(which(is.na(birthwt$low_numeric)))
}

# Create a new variable ftv2
# Collapse categories: if ftv >= 2, set to '2+', otherwise keep original value
birthwt$ftv2 <- factor(ifelse(birthwt$ftv >= 2, "2+", as.character(birthwt$ftv)))

# Check the distribution of the new variable ftv2
print(table(birthwt$ftv2))

## 
##   0   1  2+ 
## 100  47  42

# Summarize low birthweight probabilities by levels of ftv2
if (nrow(birthwt) > 0 && all(!is.na(birthwt$low_numeric))) {
  low_probabilities <- aggregate(low_numeric ~ ftv2, data = birthwt, FUN = mean)
  print(low_probabilities)
} else {
  cat("No valid rows available for aggregation.\n")
}

##   ftv2 low_numeric
## 1    0   0.3600000
## 2    1   0.2340426
## 3   2+   0.2857143

Creation of ftv2:
- Categories with ftv >= 2 are collapsed into “2+”, while smaller values (0 and 1) are retained individually.
- This simplifies analysis while addressing small sample size issues.
Summarizing Probabilities:
- The mean probability of low birthweight (low) is calculated for each level of ftv2.
- This provides insights into how physician visits during the first trimester relate to low birthweight.
  \vspace{20pt}

Using the newly created variables in d) and e), reassess the logistic regression model arrived at in b) (use and in the modeling). Comment on what you find - are the new versions of these variables important in predicting low birthweight??

# Convert categorical variables to factors with meaningful labels
birthwt$race <- factor(birthwt$race, labels = c("White", "Black", "Other"))
birthwt$smoke <- factor(birthwt$smoke, labels = c("Non-Smoker", "Smoker"))
birthwt$ht <- factor(birthwt$ht, labels = c("No", "Yes"))
birthwt$ui <- factor(birthwt$ui, labels = c("No", "Yes"))
birthwt$low <- factor(birthwt$low, labels = c("No", "Yes"))

# Create new collapsed variables ptl2 and ftv2
birthwt$ptl2 <- factor(ifelse(birthwt$ptl > 0, "1+", "0"))
birthwt$ftv2 <- factor(ifelse(birthwt$ftv >= 2, "2+", as.character(birthwt$ftv)))

# Fit logistic regression model using new variables ptl2 and ftv2
logit_model_new <- glm(low ~ age + lwt + race + smoke + ht + ui + ptl2 + ftv2,
                       data = birthwt, family = binomial())

# Perform stepwise selection to refine the model
logit_model_new_final <- step(logit_model_new, direction = "both")

## Start:  AIC=217.48
## low ~ age + lwt + race + smoke + ht + ui + ptl2 + ftv2
## 
##         Df Deviance    AIC
## - ftv2   2   196.83 214.83
## - age    1   196.42 216.42
## <none>       195.48 217.48
## - ui     1   197.59 217.59
## - smoke  1   198.67 218.67
## - race   2   201.23 219.23
## - lwt    1   200.95 220.95
## - ht     1   202.93 222.93
## - ptl2   1   203.58 223.58
## 
## Step:  AIC=214.83
## low ~ age + lwt + race + smoke + ht + ui + ptl2
## 
##         Df Deviance    AIC
## - age    1   197.85 213.85
## <none>       196.83 214.83
## - ui     1   199.15 215.15
## - race   2   203.24 217.24
## - smoke  1   201.25 217.25
## + ftv2   2   195.48 217.48
## - lwt    1   201.83 217.83
## - ptl2   1   203.95 219.95
## - ht     1   204.01 220.01
## 
## Step:  AIC=213.85
## low ~ lwt + race + smoke + ht + ui + ptl2
## 
##         Df Deviance    AIC
## <none>       197.85 213.85
## - ui     1   200.48 214.48
## + age    1   196.83 214.83
## + ftv2   2   196.42 216.42
## - smoke  1   202.57 216.57
## - race   2   205.47 217.47
## - lwt    1   203.82 217.82
## - ptl2   1   204.22 218.22
## - ht     1   205.16 219.16

# Summarize the final model
summary(logit_model_new_final)

## 
## Call:
## glm(formula = low ~ lwt + race + smoke + ht + ui + ptl2, family = binomial(), 
##     data = birthwt)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.125326   0.967561  -0.130  0.89694   
## lwt         -0.015918   0.006954  -2.289  0.02207 * 
## raceBlack    1.300856   0.528484   2.461  0.01384 * 
## raceOther    0.854414   0.440907   1.938  0.05264 . 
## smokeSmoker  0.866582   0.404469   2.143  0.03215 * 
## htYes        1.866895   0.707373   2.639  0.00831 **
## uiYes        0.750649   0.458815   1.636  0.10183   
## ptl21+       1.128857   0.450388   2.506  0.01220 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 197.85  on 181  degrees of freedom
## AIC: 213.85
## 
## Number of Fisher Scoring iterations: 4

# Obtain odds ratios for interpretation
exp_coef <- exp(coef(logit_model_new_final))
print(exp(coef(logit_model_new_final)))

## (Intercept)         lwt   raceBlack   raceOther smokeSmoker       htYes 
##   0.8822092   0.9842076   3.6724379   2.3499973   2.3787659   6.4681832 
##       uiYes      ptl21+ 
##   2.1183740   3.0921190

# Explanation:
# - We fit a logistic regression model including the newly created collapsed variables (ptl2 and ftv2).
# - Stepwise selection (step()) is used to choose important predictors based on AIC.
# - Odds ratios (exp(coef())) quantify how each predictor affects odds of low birthweight.

# Evaluate significance of new variables (ptl2 and ftv2)
summary(logit_model_new_final)

## 
## Call:
## glm(formula = low ~ lwt + race + smoke + ht + ui + ptl2, family = binomial(), 
##     data = birthwt)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.125326   0.967561  -0.130  0.89694   
## lwt         -0.015918   0.006954  -2.289  0.02207 * 
## raceBlack    1.300856   0.528484   2.461  0.01384 * 
## raceOther    0.854414   0.440907   1.938  0.05264 . 
## smokeSmoker  0.866582   0.404469   2.143  0.03215 * 
## htYes        1.866895   0.707373   2.639  0.00831 **
## uiYes        0.750649   0.458815   1.636  0.10183   
## ptl21+       1.128857   0.450388   2.506  0.01220 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 234.67  on 188  degrees of freedom
## Residual deviance: 197.85  on 181  degrees of freedom
## AIC: 213.85
## 
## Number of Fisher Scoring iterations: 4

# Explanation:
# - Check p-values for ptl2 and ftv2 in the summary output.
# - If these variables have significant p-values (typically < 0.05), they are important predictors.

How to interpret your findings:

After running this code, inspect the summary output:
- If ptl2 or ftv2 have significant p-values (typically < 0.05), it indicates that these new collapsed versions are important predictors of low birthweight.
- If they are not significant, it suggests that collapsing categories did not substantially improve their predictive utility.
  
  \vspace{20pt}

In a manner similar to the approach used in the book, split the data into a training and test set, where the test set is about 20% the size of the entire dataset. Then, using variables that are justifiable for inclusion in discriminant analysis, fit LDA and QDA models to the training set and form confusion matrices, calculate the sensitivity, specificity, and the accuracy of each method using the test set, and do the same for the logistic regression models built in f) and b). Which model performs the best? Remember you MUST set the seed using the package in a manner similar to as done in the notes (but don’t use my name to set the seed!)

# Set seed using TeachingDemos package
char2seed("birthwt_analysis")

# Split the data into training (80%) and test (20%) sets
set.seed(123)  # Ensure reproducibility
train_indices <- sample(seq_len(nrow(birthwt)), size = 0.8 * nrow(birthwt))
train_data <- birthwt[train_indices, ]
test_data <- birthwt[-train_indices, ]

# Check the dimensions of the training and test sets
cat("Training set dimensions:", dim(train_data), "\n")

## Training set dimensions: 151 13

cat("Test set dimensions:", dim(test_data), "\n")

## Test set dimensions: 38 13

# Fit LDA model on training data
lda_model <- lda(low ~ age + lwt + race + smoke + ht + ui + ptl2 + ftv2, data = train_data)
lda_pred <- predict(lda_model, newdata = test_data)$class

# Fit QDA model on training data
qda_model <- qda(low ~ age + lwt + race + smoke + ht + ui + ptl2 + ftv2, data = train_data)
qda_pred <- predict(qda_model, newdata = test_data)$class

# Fit Logistic Regression model using new variables (from part f)
logit_model <- glm(low ~ age + lwt + race + smoke + ht + ui + ptl2 + ftv2, data = train_data, family = binomial())
logit_prob <- predict(logit_model, newdata = test_data, type = "response")
logit_pred <- ifelse(logit_prob > 0.5, "Yes", "No")

# Evaluate performance for each model
actual <- test_data$low

# Confusion matrices
lda_cm <- table(Predicted = lda_pred, Actual = actual)
qda_cm <- table(Predicted = qda_pred, Actual = actual)
logit_cm <- table(Predicted = logit_pred, Actual = actual)

cat("Confusion Matrix for LDA:\n")

## Confusion Matrix for LDA:

print(lda_cm)

##          Actual
## Predicted No Yes
##       No  23   9
##       Yes  4   2

cat("\nConfusion Matrix for QDA:\n")

## 
## Confusion Matrix for QDA:

print(qda_cm)

##          Actual
## Predicted No Yes
##       No  22   7
##       Yes  5   4

cat("\nConfusion Matrix for Logistic Regression:\n")

## 
## Confusion Matrix for Logistic Regression:

print(logit_cm)

##          Actual
## Predicted No Yes
##       No  24   9
##       Yes  3   2

# Calculate sensitivity, specificity, and accuracy for each model
calc_metrics <- function(cm) {
  sensitivity <- cm["Yes", "Yes"] / sum(cm[, "Yes"])
  specificity <- cm["No", "No"] / sum(cm[, "No"])
  accuracy <- sum(diag(cm)) / sum(cm)
  return(list(Sensitivity = sensitivity, Specificity = specificity, Accuracy = accuracy))
}

lda_metrics <- calc_metrics(lda_cm)
qda_metrics <- calc_metrics(qda_cm)
logit_metrics <- calc_metrics(logit_cm)

cat("\nPerformance Metrics for LDA:\n")

## 
## Performance Metrics for LDA:

print(lda_metrics)

## $Sensitivity
## [1] 0.1818182
## 
## $Specificity
## [1] 0.8518519
## 
## $Accuracy
## [1] 0.6578947

cat("\nPerformance Metrics for QDA:\n")

## 
## Performance Metrics for QDA:

print(qda_metrics)

## $Sensitivity
## [1] 0.3636364
## 
## $Specificity
## [1] 0.8148148
## 
## $Accuracy
## [1] 0.6842105

cat("\nPerformance Metrics for Logistic Regression:\n")

## 
## Performance Metrics for Logistic Regression:

print(logit_metrics)

## $Sensitivity
## [1] 0.1818182
## 
## $Specificity
## [1] 0.8888889
## 
## $Accuracy
## [1] 0.6842105

# Compare models based on metrics
best_model <- which.max(c(lda_metrics$Accuracy, qda_metrics$Accuracy, logit_metrics$Accuracy))
model_names <- c("LDA", "QDA", "Logistic Regression")
cat("\nBest Model Based on Accuracy:", model_names[best_model], "\n")

## 
## Best Model Based on Accuracy: QDA

Explanation:

Data Splitting:
- The dataset is split into training (80%) and test (20%) sets using sample().
- The seed is set using char2seed() from the TeachingDemos package to ensure reproducibility.
Model Fitting:
- LDA (lda()), QDA (qda()), and logistic regression (glm()) models are fit using the training data.
- Predictions are made on the test set.
Confusion Matrices:
- Confusion matrices are created for each model to compare predicted vs actual values.
Performance Metrics:
- Sensitivity: Proportion of correctly predicted positive cases.
- Specificity: Proportion of correctly predicted negative cases.
- Accuracy: Overall proportion of correct predictions.
Comparison:
- Metrics are calculated for all models.
- The model with the highest accuracy is identified as the best-performing model.

Output:

Confusion matrices for LDA, QDA, and logistic regression.
Sensitivity, specificity, and accuracy for each model.
Identification of the best-performing model based on accuracy.
\vspace{20pt}

Using your final model from f), interpret the estimates for all covariates.

# Convert categorical variables to factors with meaningful labels
birthwt$race <- factor(birthwt$race, labels = c("White", "Black", "Other"))
birthwt$smoke <- factor(birthwt$smoke, labels = c("Non-Smoker", "Smoker"))
birthwt$ht <- factor(birthwt$ht, labels = c("No", "Yes"))
birthwt$ui <- factor(birthwt$ui, labels = c("No", "Yes"))
birthwt$low <- factor(birthwt$low, labels = c("No", "Yes"))

# Create new collapsed variables ptl2 and ftv2
birthwt$ptl2 <- factor(ifelse(birthwt$ptl > 0, "1+", "0"))
birthwt$ftv2 <- factor(ifelse(birthwt$ftv >= 2, "2+", as.character(birthwt$ftv)))

# Fit the final logistic regression model using new variables (from part f)
logit_model_final <- glm(low ~ age + lwt + race + smoke + ht + ui + ptl2 + ftv2,
                         data = birthwt, family = binomial())

# Get the coefficients from the model
coefficients <- coef(logit_model_final)

# Exponentiate the coefficients to calculate odds ratios
odds_ratios <- exp(coefficients)

# Display coefficients and odds ratios in a table
results <- data.frame(
  Variable = names(coefficients),
  Coefficient = coefficients,
  Odds_Ratio = odds_ratios
)
print(results)

##                Variable Coefficient Odds_Ratio
## (Intercept) (Intercept)  0.82301896  2.2773647
## age                 age -0.03723429  0.9634504
## lwt                 lwt -0.01565301  0.9844689
## raceBlack     raceBlack  1.19241321  3.2950232
## raceOther     raceOther  0.74068488  2.0973715
## smokeSmoker smokeSmoker  0.75552837  2.1287360
## htYes             htYes  1.91316586  6.7745020
## uiYes             uiYes  0.68019548  1.9742636
## ptl21+           ptl21+  1.34376338  3.8334431
## ftv21             ftv21 -0.43637967  0.6463723
## ftv22+           ftv22+  0.17900852  1.1960309

# Explanation:
# - The `coef()` function extracts the coefficients (log-odds) from the logistic regression model.
# - The `exp()` function transforms these coefficients into odds ratios.
# - Odds ratios represent the multiplicative change in odds of low birthweight for a one-unit increase in a predictor (or a change from reference level for categorical predictors).

# Interpretation of results:
cat("\nInterpretation of Results:\n")

## 
## Interpretation of Results:

for (i in seq_along(results$Variable)) {
  var <- results$Variable[i]
  coef_val <- results$Coefficient[i]
  or_val <- results$Odds_Ratio[i]
  
  if (var == "(Intercept)") {
    cat(sprintf("- The intercept represents the baseline log-odds when all predictors are at their reference levels. Odds ratio: %.3f\n", or_val))
  } else {
    if (coef_val > 0) {
      cat(sprintf("- A one-unit increase in '%s' increases the odds of low birthweight by a factor of %.3f.\n", var, or_val))
    } else if (coef_val < 0) {
      cat(sprintf("- A one-unit increase in '%s' decreases the odds of low birthweight by a factor of %.3f.\n", var, or_val))
    } else {
      cat(sprintf("- '%s' has no effect on the odds of low birthweight.\n", var))
    }
  }
}

## - The intercept represents the baseline log-odds when all predictors are at their reference levels. Odds ratio: 2.277
## - A one-unit increase in 'age' decreases the odds of low birthweight by a factor of 0.963.
## - A one-unit increase in 'lwt' decreases the odds of low birthweight by a factor of 0.984.
## - A one-unit increase in 'raceBlack' increases the odds of low birthweight by a factor of 3.295.
## - A one-unit increase in 'raceOther' increases the odds of low birthweight by a factor of 2.097.
## - A one-unit increase in 'smokeSmoker' increases the odds of low birthweight by a factor of 2.129.
## - A one-unit increase in 'htYes' increases the odds of low birthweight by a factor of 6.775.
## - A one-unit increase in 'uiYes' increases the odds of low birthweight by a factor of 1.974.
## - A one-unit increase in 'ptl21+' increases the odds of low birthweight by a factor of 3.833.
## - A one-unit increase in 'ftv21' decreases the odds of low birthweight by a factor of 0.646.
## - A one-unit increase in 'ftv22+' increases the odds of low birthweight by a factor of 1.196.

Fit Logistic Regression Model:
- The final logistic regression model includes predictors such as age, lwt, race, smoke, ht, ui, and the newly created variables ptl2 and ftv2.
Extract Coefficients:
- The coefficients from the model represent changes in log-odds for a one-unit increase in continuous predictors or a change from the reference level for categorical predictors.
Calculate Odds Ratios:
- Exponentiating the coefficients (exp()) converts them into odds ratios, which are easier to interpret.
Interpretation:
- For each predictor:
  - A positive coefficient means an increase in odds of low birthweight.
  - A negative coefficient means a decrease in odds.
  - An odds ratio greater than 1 indicates increased odds; less than 1 indicates decreased odds.

Akshay Kumar Exam1 (Take Home)

INFO-H 415/515

03/12/2025

Please hand this result in on Canvas no later than 11:59pm on Wednesday, March 12! Do not work in groups!

Linear Regression with cpus Data

Output:

Key Observations:

Key Points:

How to interpret your findings:

Explanation:

Output: