options(repos = "https://cran.rstudio.com/")
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(caret)
## Loading required package: ggplot2
## Loading required package: lattice
library(ModelMetrics)
## 
## Attaching package: 'ModelMetrics'
## The following objects are masked from 'package:caret':
## 
##     confusionMatrix, precision, recall, sensitivity, specificity
## The following object is masked from 'package:base':
## 
##     kappa
library(datarium)
library(corrplot)
## corrplot 0.92 loaded
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(bestglm)
## Loading required package: leaps
library(glmnet)
## Loading required package: Matrix
## Loaded glmnet 4.1-8

Part 1

Sample n rows from a table

data("marketing", package = "datarium")
sample_n(marketing, 3)  
##   youtube facebook newspaper sales
## 1  252.84    35.40     11.16 22.08
## 2  104.64    14.16     31.08 12.72
## 3  226.08    21.72     30.72 17.88

##Sales Distribution

p <- ggplot(marketing) +
  geom_histogram(aes(x = sales, y = ..density..),
                 binwidth = 1, fill = "grey", color = "black") + geom_density(aes(x=sales, color="red"),
                 show.legend = FALSE)
p + theme_bw()
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Scaling Techniques (numeric variables)

Standardization \(z_i=\frac{x_i-\bar{x}}{\sigma}\)

preproc1 <- preProcess(marketing, method=c("center", "scale"))
norm1 <- predict(preproc1, marketing)
head(norm1, 5)
##       youtube   facebook newspaper      sales
## 1  0.96742460  0.9790656 1.7744925  1.5481681
## 2 -1.19437904  1.0800974 0.6679027 -0.6943038
## 3 -1.51235985  1.5246374 1.7790842 -0.9051345
## 4  0.05191939  1.2148065 1.2831850  0.8581768
## 5  0.39319551 -0.8395070 1.2785934 -0.2151431

Min-Max Scaling \((x-x_{min}/({x_{max}-x_{min})}\)

preproc2 <- preProcess(marketing, method=c("range"))
norm2 <- predict(preproc2, marketing)
head(norm2, 5)
##     youtube  facebook newspaper     sales
## 1 0.7757863 0.7620968 0.6059807 0.8070866
## 2 0.1481231 0.7923387 0.3940193 0.3464567
## 3 0.0557998 0.9254032 0.6068602 0.3031496
## 4 0.5099763 0.8326613 0.5118734 0.6653543
## 5 0.6090632 0.2177419 0.5109938 0.4448819
M <-cor(norm1)
p.mat <- cor.mtest(norm1)
corrplot(M, type="upper", order="hclust",
         p.mat = p.mat$p, sig.level = 0.05)

library caret - split in 80% training data with createDataPartition(), where y = a vector of outcomes p = the percentage of data that goes to training list = logical should the results be in a list (TRUE) or a matrix (FALSE)

set.seed(123)
training.samples <- createDataPartition(y = norm1$sales, p = 0.8, list = FALSE)
train.data  <- norm1[training.samples, ]
test.data <- norm1[-training.samples, ]
model <- lm(sales ~., data = train.data)
predictions <- predict(model,test.data)

Model Performance

data.frame( RMSE = RMSE(predictions, test.data$sales),
R2 = R2(predictions, test.data$sales),
MAE = MAE(predictions, test.data$sales),
MSE = mse(predictions, test.data$sales))
##        RMSE        R2       MAE        MSE
## 1 0.3139314 0.9049049 0.2289764 0.09855291

multicollinearity VIF if VIFs > 4 - a cause for concern if VIFs = 10 or above - a very high degree of multicollinearity If VIF is high, remove that variable

vif(model)
##   youtube  facebook newspaper 
##  1.004440  1.118155  1.115449

Part 2

Linear Regression

model_swiss = lm(Fertility ~ .,data = swiss)
lm_coeff = model_swiss$coefficients
print(summary(model_swiss))
## 
## Call:
## lm(formula = Fertility ~ ., data = swiss)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.2743  -5.2617   0.5032   4.1198  15.3213 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
## Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
## Examination      -0.25801    0.25388  -1.016  0.31546    
## Education        -0.87094    0.18303  -4.758 2.43e-05 ***
## Catholic          0.10412    0.03526   2.953  0.00519 ** 
## Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.165 on 41 degrees of freedom
## Multiple R-squared:  0.7067, Adjusted R-squared:  0.671 
## F-statistic: 19.76 on 5 and 41 DF,  p-value: 5.594e-10

Logistic Regression

data("SAheart")
model_cholecterol = glm(chd ~ ldl, data = SAheart, family = binomial)

plot(jitter(chd, factor = 0.1) ~ ldl, data = SAheart, pch = 20, 
     ylab = "Probability of CHD", xlab = "Low Density Lipoprotein Cholesterol")
grid()

curve(predict(model_cholecterol, data.frame(ldl = x), type = "response"), 
      add = TRUE, col = "dodgerblue", lty = 2)

summary(model_cholecterol)
## 
## Call:
## glm(formula = chd ~ ldl, family = binomial, data = SAheart)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.96867    0.27308  -7.209 5.63e-13 ***
## ldl          0.27466    0.05164   5.319 1.04e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 596.11  on 461  degrees of freedom
## Residual deviance: 564.28  on 460  degrees of freedom
## AIC: 568.28
## 
## Number of Fisher Scoring iterations: 4

Ridge Regression

Regularization is generally useful in the following situations: - Large number of variables - Low ratio of number observations to number of variables - High Multi-Collinearity

X = swiss[,-1]
y = swiss[,1]
  • Use cv.glmnet( ) function to do cross validation
  • Use alpha = 0 for ridge regression
  • Use the best lambda by using lambda.min
set.seed(123) 

model_ridge = cv.glmnet(as.matrix(X),y,alpha = 0,lambda = 10^seq(4,-1,-0.1))

best_lambda = model_ridge$lambda.min
ridge_coeff = predict(model_ridge,s = best_lambda,type = "coefficients")
ridge_coeff
## 6 x 1 sparse Matrix of class "dgCMatrix"
##                           s1
## (Intercept)      62.97585936
## Agriculture      -0.09863022
## Examination      -0.33967990
## Education        -0.64733678
## Catholic          0.07703325
## Infant.Mortality  1.08821833

Lasso Regression

Lasso stands for Least Absolute Shrinkage and Selection Operator. - Use the same swiss dataset and X and Y - Use glmnet for cross-validation - Set standartize = TRUE (this is default)

set.seed(123)
model = cv.glmnet(as.matrix(X),y,alpha = 1,lambda = 10^seq(4,-1,-0.1))
best_lambda = model$lambda.min
lasso_coeff = predict(model,s = best_lambda,type = "coefficients")
lasso_coeff 
## 6 x 1 sparse Matrix of class "dgCMatrix"
##                           s1
## (Intercept)      65.46374579
## Agriculture      -0.14994107
## Examination      -0.24310141
## Education        -0.83632674
## Catholic          0.09913931
## Infant.Mortality  1.07238898

Both ridge regression and lasso regression are addressed to deal with multicollinearity.