R notes

Question 1

library(dplyr)

Answer 1

pipes, group by, summarise

group_by(factor column)

summarise( Means = mean(), Medians = median(), Freq = n())

Question 1

Combining levels

Answer 1

library(plyr)
mapvalues(dataset$column, levels(dataset$column), c(“name level1”, etc.)

Question 2

Unique values

Answer 2

length(unique(dataset$column))

Question 3

Create Factor of intervals

Answer 3

cut(dataset$column,
breaks = c(0, 15, 35, 65, 100)
labels = c(“Child”, “Young Adult”, Middle.aged”, “Elderly”))

Question 4

Binary Target - Factor to Numeric

Answer 4

ifelse(factor == “success level”, 1, 0)

Question 5

Relevel

Answer 5

levels(factor)

relevel(dataset$column, ref = “new reference level”)

Question 6

Scatterplot variation

Answer 6

ggplot(dataset, aes(x = variable, y = target)) +
geom_count()

Question 7

Regularization- code

Answer 7

requires matrix object for x: Xmat <- model.matrix(target ~ predictors, data = dataset)
library(glmnet)
rid.mod1 <- glmnet( x = Xmat[, -1], y = dataset$target,
family = “gaussian”,
lambda = ?,
alpha = ?,
standardize = T <- default, so not needed
)

coef(rid.mod1) #all estimates in one place
plot(cv.glmet()) - produces a CV visual; the left dashed line indicates where CV error is minimized; the right dashed line indicates the one-standard-error rule; the numbers at the top indicate the number of predictors in the model

Question 8

Regularization - notes

Answer 8

hyperparameters: alpha and lambda
-alpha = 0 -> Ridge; alpha = 1 -> Lasso
-lambda = 0 -> no shrinkage

predict(rid.mod4, newx = Xmat[, -1])
training RMSE <- sqrt(sum((dataset$target - predict(rid.mod4, newx = Xmat[, -1]))^2) / nrow(dataset)

Question 9

Symbols

Answer 9

~ = ‘explained by’
+ = ‘include’
. = ‘all remaining variables’
- = ‘exclude’
: = ‘interact’
* = ‘cross’ -> include and interact

Question 10

Regularization - CV code

Answer 10

set.seed(161)
rid.cv <- cv.glmnet( x = Xmat[, -1], y = dataset$target,
family = “gaussian”,
alpha = ?,
lambda = ??,
standardize = T <- default, so not needed
)
-autogenerates lambda if not specified
-nfold is 10 by default

rid.cv$lambda[, ] # Shows results for the 100 lambda values (100 is default)
rid.cv$lambda.min #returns best lambda OR rid.cv$lambda.1se
–>use rid.cv$lambda.min for lambda in regular glmnet model

coef(rid.mod1)

Question 11

Train function pt 1

Answer 11

library(caret)
set.seed(161)
tree.cv1 <- train(
target ~ predictors,
data = dataset,
method = “rpart”,
trControl = trainControl(method = “cv”, number = 10),
metric = “RMSE”,
tuneGrid = expand.grid(cp = seq(0.0, 0.005, 0.0005)),
control = rpart.control(minsplit = , minbucket = maxdepth = ),
na.action = na.pass
)

plot(tree.cv1)
tree.cv1$results
tree.mod# <- tree.cv1$finalModel #extract tree with lowest cv
rpart.plot(tree.mod#)
tree.cv1$bestTune #extract the optimal cp value

control optional, will use defaults
na.action shouldn’t be needed if we adequately handle missing data at the start
can replace 1st 2 arguments with:
y = dataset$column, x = dataset[, c(“predictor1”, “predictor2”, etc.)]

Question 12

Train pt 2

Answer 12

need a binary factor target
library(caret)
set.seed(161)
tree.cv1 <- train(
factor(target) ~ predictors,
data = dataset,
method = “rpart”,
trControl = trainControl(method = “cv”, number = 10),
metric = “Accuracy”,
tuneGrid = expand.grid(cp = seq(0.0, 1, 0.01)),
na.action = na.pass,
parms = list(split = “gini” or “information”)
)

plot(tree.cv1) #higher the accuracy the better

tree.mod# <- tree.cv1$finalModel #extract the model with the highest cv accuracy
rpart.plot(tree.mod#, extra = 4) #makes tree
tree.cv$bestTune #extract optimal cp value

“Accuracy” - %age of correct predictions = 1 - classification error rate

Question 13

Other

Answer 13

poly(variable, 2) = variable + I(variable^2)

Question 14

Train pt 3

Question 15

RF plots

Answer 15

library(pdp)

partial(rf.mod#, pred.var = “variable”, plot = T)

Question 16

PCA

Answer 16

base package
pca1 <- promp(~. - factor, data = dataset, center = T, scale. = T)
OR pca1 <- prcomp(dataset[, colnames(dataset) != “factor”], center = T, scale. = T)
-exclude factors because this function doe not work with factors
-center default = T

pca1$rotation #access the loadings
pca1$x #obtain scores of the PCs (for each observation)
summary(pca1) # proportion of variance explained
plot(pca1, type = “line”) #scree plot
OR screeplot(pca1, type = “line”, npcs = # of PCs to plot)

pca1$x[, 1 or “PC1”] #scores of PC1
dataset.new <- dataset
dataset.new$PC1 <- pca1$x[, “PC1”]
dataset.new$variable1, 2, etc. <- NULL

Question 17

k-means clustering

Answer 17

set.seed(161)
dataset.scaled <- data.frame(scale(iris[, columns(dataset) != “factor”]))
km1 <-kmeans(
x = dataset[, c(“variable1”, “varaible2”], # Required - must be scaled
centers = 5, #Required - k = the number of clusters
iter.max = , # max iterations with default = 10, can try higher (30, etc.)
nstart = # # of random initial cluster assignments with default = 1 -> ‘best’ chosen
algorithm , #irrelevant, we viewed Lloyd which is not the default
)

km1$cluster # assignments for each obs
km1$tot.withinss #extract total within-cluster variation

km1$betweenss/dm1$totss

library(ggplot2)
ggplot(dataset.scaled, aes(x = variable1, y = variable2, color = factor(km1$cluster))) +
geom_point()

Question 18

hierarchical clustering

Answer 18

nc1 <- hclust(

plot(hc1) #dendrogram
cuttree(hc1, h = 4) # cut at height of 4
OR currtree(hc1, k = 3) # cut into 3clusters

Question 19

R code

Answer 19

lm = linear models
mlr.mod1 <- lm(target ~ predictors, data = dataset)
summary(mlr.mod1)

predict(mlr.mod1) -> y^’s
predict(mlr.mod1, newdata = )
->newdata does require specific columns but the model only uses ones included in the model formula

RMSE: sqrt(sum(dataset$target - predict(mlr.mod1, newdata = dataset)^2) / nrow(dataset))

Question 20

R code - Binarization

Answer 20

bin.mod <- dummyVars(~factor, data = dataset, fullRank = T)
predict(bin.mod, newdata = dataset)

new.dataset <- cbind(dataset[, colnames(dataset) != “factor”], predict(bin.mod, newdata = dataset))

Question 21

R Code - Stratified Sampling

Answer 21

R:
library(caret)
set.seed(161)
ind <- createDataPartition(dataset$target, p = 0.7, list = F)
train <- dataset[ind, ]
test <- dataset[-ind, ]

rm(ind) #remove… ->not necessary

mean(train$target)
mean(test$target)

Question 22

R Code: Stepwise Selection

Answer 22

Null model: mlr.mod0 <- lm(target ~ 1, data = dataset)
Biggest model: mlr.mod4 <- lm(target ~ ., data = dataset) <-largest model you want to consider

library(MASS)
AIC:
stepAIC(Null model, direction = “forward”, scope = list(upper = largest model), k = 2)

AIC -> k = 2; BIC -> k = log(nrow(dataset))
Backward–>no scope needed, and if scope is not included, the default for direction is “backward”

k = 2 (default)
direction = “backward” (default)

drop1() -> what does 1st round with backward selection and AIC produce for AIC if drop each predictor

Question 23

R code

Answer 23

function: glm

glm.mod1 <- glm(target ~ predictors, family = “distribution”(link = “link”), data = dataset, weights = column in dataset of weights)

->if using weights, target must be on a per weight basis; also exclude column of weights from predictors; only with normal, gamma, or inverse gaussian on exam

predict(glm.mod1, type = “response”) -> this returns mu^ not (1/mu^), etc.

Question 24

R code

Answer 24

use family = “binomial” and link = “logit”

exp(coef(modelname))

predict(modelname, type = “response”) -predicted probabilities NOT predicted odds

To get target predictions: ifelse(predict(modelname, type = “response”)>0.5,1,0)
-0.5 is cutoff that we select to separate the target predictions

Question 25

R code - cont.

Answer 25

library(pROC)

target <- target vector
probs <- probs vector
preds <- preds vector

roc <- roc(target, probs)
auc(roc) - returns AUC of the roc object
plot(roc) - graph of roc curve

library(caret)
confusionMatrix(factor(preds), factor(target), positive = “1”)

Question 26

R code

Answer 26

poi.mod1 <- glm(target ~ predictors, family = poisson(link = “log”), data = dataset)
exp(coeff(poi.mod1))

poi.mod2 <- glm(target ~ predictors - w, family = poisson(link = “log”), offset = log(w), data = dataset)
->w is exposure column

predict(poi.mod2, type = “response”)

Question 27

R code

Answer 27

rpart() - key hyperparameters:
-minsplit - minimum number of observations in a node
-minbucket - minimum number of observations permitted for a terminal node
-maxdepth - maximum depth of terminal nodes
-cp - complexity parameter (0 most flexible; 1 least flexible)
-xval = # of cross validation to perform

anova is default for method
maxcompete - most # of competitor splits stored for each split
maxsurrogate - most # of surrogate splits stored for each split; 0 increases speed
usesurrogate - 0 - not used; 1 - used in sequence and left behind if fails; 2 used in sequence, majority if still fails (defualt)

library(rpart)
set.seed(161)
tree.mod1 <- rpart(target ~ predictors, data = dataset, method = “anova” or “poisson”, control = rpart.control(minsplit = #, minbucket = #, maxdepth = #, cp = ,
xval = , maxcompete = , maxsurrogate = , usesurrogate = )
tree.mod1 # provides data on all nodes in the tree
summary(tree.mod1)

library(rpart.plot)
rpart.plot(tree.mod1, digits = #) default digits = 2
predict(tree.mod1)

Poisson
predict(tree.mod1) -> predicted poisson rates, not counts
predict(tree.mod1)*train$w -> obtain predicted poisson counts

Question 28

R code - pruning

Answer 28

rpart already performs cv

tree.mod1$cptable #returns cp table
columns:
cp
nsplits
rel error (relative error) -> (SSE w n splits/ SSE w 0 splits)
xerror (cv error)
x std (cv std)
CP decreases as nsplits increases
As nsplits increase, relative error decreases

cp.min <- tree.mod1$cptable[whichmin(tree.mod1$cptable[, “xerror]),”CP”)

tree.mod1a <- prune.rpart(tree.mod1, cp = cp.min

plotcp(tree.mod1) # can select cp with leftmost point below the horizontal dotted line (one standard error rule)

Question 29

R code

Answer 29

need binary target
library(rpart)
set.seed(161)
tree.mod# <- rpart(target ~ predictors, data = dataset, method = “class”, parms = list(split = “gini” or “information”))
library(rpart.plot)
rpart.plot(tree.mod#, extra = 4) #4 displays probs for each class @ each node
predict(tree.mod#, type = “prob” or “class”)
#prob - explicitly display probs (default), class - explicitly display target predictions

Extract only the positive class probabilities - for use in ROC curve
predict(tree.mod#, type = “prob)[, 2] or [, “1”]

tree.mod#
node), split, n, loss (# misclassified observations), yval (prediction), (yprob) (probabilities for each class)

Question 30

R code - pruning

Answer 30

tree.mod#$cptable #complexity parameter table
-cp
-nsplit
-relerror (classification error rate with n split/classification error w 0 splits)
-xerror (cverror)
-xstd

classification error rate = %age of observations with wrong predictions
If cv error tied, we favor less flexible model due to parsimony
One-standard-error rule -> leftmost point below the horizontal line

plotcp(tree.mod#)
cp.min <- tree.mod#$cptable[which.min(tree.mod#$cptable[, “xerror”], “CP”]
tree.mod#+1 <- prune.rpart(tree.mod#, cp = cp.min)
rpart.plot(tree.mod#+1, extra = 4)

Question 31

R code - Random Forests

Answer 31

library(randomForest)
set.seed(161)
rf.mod1 <- randomForest(
target ~ predictors,
data = dataset,
importance = T, # whether to examine the importance of features
ntree = 101, # number of trees
mtry = 2, # number of candidate predictors for each split
nodesize = 5, # min number of observations permitted for a terminal node
maxnodes = #, # maximum number of terminal nodes
keep.forest = T # if want to obtain predictions for all trees later
)

rf.mod1$importance #variable importance measure
OR
library(caret)
varImp(rf.mod1)

predict(rf.mod1) #predictions on out-of-bag obs
->for all trees where that observation is out of bag and averages. If one is NOT out-of-bag for any tree, then NA
predict(rf.mod1, newdata = dataset) # to get predictions for all obs
predict(rf.mod1, newdata = dataset, predict.all = T)$individual/aggregate

mtry default: Integer and >=1
Factor -> floor(sqrt(ncol(x))) sqrt(predictors)
Not Factor -> max(floor(ncol(x)/3),1) predictors / 3

Question 32

R code - Boosting (gbm)

Answer 32

library(gbm)
gbm.mod1 <- gbm(
target ~ predictors,
data = dataset,
distribution = “guassian/bernoulli/poisson”, # R will guess if empty
n.trees = 300, # number of trees
interaction.depth = 2, # maximum depth of terminal nodes
shrinkage = 0.01, # shrinkage parameter
n.minobsinnode = 10 # min number of obs permitted for a terminal node
bag.fraction = 1, # portion of observations used for each tree if <1, set seed
)

summary(gmb.mod1) #variable importance measure
predict(gbm.mod1, n.trees = 1) #default n.trees = all

Question 33

R code - Boosting (xgboost)

Answer 33

library(xgboost)
requires new data type for datasets
m.f <- model.frame(target ~ predictors, data = head(dataset)
m.m <- model.matrix(attr(m.f., “terms”), data = dataset)
d.m <- xgb.DMatrix(m.m, label = dataset$target)

xgb.parm1 <- list(
booster = “gbtree”, #decision trees
objective = “reg:squarederror/binary:logistic”, #regression/classification
eval.metric = “rmse/auc”, #regression/classification
max.depth = 2,
eta = 0.01, # shrinkage parameter
gamma = 0, #min reduction of splitting measure
subsample = 1, #proportion of data used
colsample_bytree = 1, #proportion of predictors used
min_child_weight = 10, #min obs in node for regression; purity of node for class.
)

xgb.mod1 <- xgb.train(data = d.m, nrounds = 300, params = xgb.parm1)
nrounds = xgb.cv1$bestiteration #extracts best nrounds

set.seed(161)
xgb.cv1 <- xgb.cv(data = d.m, nrounds = 2000, params = xgb.parm1,
nfold = 5,
prediction = F,
print_every_n = 100,
early_stopping_rounds = 50,
maximize = F
)
#last line is optimal value for nrounds

xgb.importance(model = xgb.mod1)
#column gain = %age of contribution to the boosted model

predict(xgb.mod1, newdata = d.m)

Question 34

R code

Answer 34

OVER
in.over <- c(which(train$y == 0), rep(which(train$y == 1), times = k))
train.over <- train[ind.over, ]

UNDER
set.seed(161)
ind.under <- c(which(train$y == 1), sample(which(train$y == 0), size = sum(train$y == 1)))
train.under <- train[ind.under, ]

Question 35

Bivariate - numerical in R

Answer 35

cor(dataset$column, dataset$column) or cor(dataset[,~factors])
table(dataset$column, dataset$column)

library(dplyr)
dataset %>%
group_by(factor) %>%
summarise(Means = mean(target),
Medians = median(target),
Freq = n())

Question 36

Univariate - graphical in R

Answer 36

library(ggplot2)
ggplot(data = dataset, mapping = aes(x = variable)) +
* geom_histogram/bar/boxplot() +
labs(x = “label name”)

*
histogram - (bins = ) #bins optional but can adjust
bar - (position = “dodge”) or (position = “fill”) # fill = area of 1
boxplot - (position = “dodge”) or (position = “fill”) # fill = area of 1

Question 37

Bivariate - graphical in R (histogram and bar)

Answer 37

library(ggplot2)
ggplot(data = dataset, mapping = aes(x = variable, fill = factor)) +
* geom_histogram/bar() +
labs(x = “label name”)

*
histogram - (position = “dodge”, bins = ) #bins optional but can adjust
bar - (position = “dodge”) or (position = “fill”) # fill = area of 1
Consider density instead of fill
-> add “y = ..density..” after fill statement on line 1

Question 38

Bivariate - graphical in R (scatterplot and boxplot and facet)

Answer 38

library(ggplot2)
ggplot(data = dataset, mapping = aes(x = variable, y = target)) +
geom_point/count/boxplot() +
labs(x = “label name”, y = “label name”)

ggplot(data = dataset, mapping = aes(x = variable, y = target)) +
geom_point() +
labs(x = “label name”, y = “label name”) +
* facet_wrap(~factor)

*
scales = “free” - if don’t want same scaling