Formulario Parte MAT dell'esame di Business data analytics

#LIBRERIE

#lab1

library(ggplot2)

library(ggthemes)

library(GGally)

#lab2

library(tidyverse) #anova

#lab3

library(mvtnorm) #PCA

library(ggfortify) #PCA -> autoplot() -> scatterplot of scores #useful also for Lab7

library(rgl)

library(pracma)

#lab4

library(rgl)

library(cluster)

#lab5

library(MASS) #for stepwise forward selection too

library(car)

library(rgl)

#lab6

install.packages("boot")

library(boot) #cross-validation

#lab7

library(ISLR)

library(pls) #pca automated with CV

library(glmnet) #ridge and lasso

library(car)

#lab8

library(PRROC) #ROC Curve

library(class) #knn

#lab7

install.packages("tree")

library(tree) #tree

install.packages("randomForest")

library(randomForest)

#GES

#NETWORK ANALYSIS

library(igraph)

library(dplyr)

library(anytime)

library("reshape2")

library(RColorBrewer)

library(ggplot2)

library(corrplot) #for correlation matrix

#INFERENCE

library(MatchIt) #PSM

library(ivreg) #IV

library(MASS)

library(mvtnorm)

library(rdrobust) #SHARP RDD

library(rddensity)

library("poweRlaw")

#LABORATORIO 1# (BASIS)

install.packages(c("ggplot2", "ggthemes"))

library(ggplot2)

library(ggthemes)

rm(list=ls()) #cancella tutto enviroment

cat("\014") #pulire console

x11()

graphics.off()

rm(x) #cancella oggetto x in enviroment

getwd() #stampa cartella di lavoro / current working directory

dir() #visualizza file disponibili nella cartella di lavoro

dev.off() #close all graphs

dati_... = read.csv(".txt") # apre dataset .txt e gli dà il nome sui puntini (di base campi

separati da virgole)

dati_... = read.csv(".txt", sep =" ") #apre dataset con campi separati da spazio

head(dati_...) #stampa le prime righe del dataset

dati_...$... #per accedere ad una colonna in particolare del dataset

dati_...$... = NULL #cancello una colonna in particolare

str(dati_...) #dataframe structure: stampa tutte le covariate e mi dice cosa sono

dati_pinguini$specie = factor(dati_pinguini$specie) #es. per fattorizzare variabili

categoriche (quando ho variabili chr devo sempre farlo)

#serve per specificare a R che si tratta di variabili categoriche (che quindi possono

assumere solo certi valori finiti)

################################################

###### Creation of Numeric Vectors

###############################################

v1 = c(1, 2.4, 5.6, 7)

v1[3]

v2 = c(2,3,4,0)

## Ci sono molti altri metodi per creare vettori. Uno molto utile è il seguente

v3 = seq(1,10, by = 0.1) # Numeri equispaziati di 0.1, da 1 a 10

v3

v3.bis = seq(1, 10, length.out=9) #Numbers from 1 to 10, equally spaced so that

#the resulting vector has 9 elements

v3.bis

## Principal Vectors' Operations

2 * v1 # Usual product by a number

v1/2 # Usual division by a number

v1 - 1 #1 is subtracted from all elements

v1 + 4 # 4 is addedo to all elements

v1 + v2 # Sum Component by Componente

v1 * v2 # Product componeny by Component

sum(v1) # Sum of all elements of V1

v1[2] # Second element of V1

v1[c(2,4)] #Second and fourth element of v1

v1[1:3] # Elements of v1 from first to third

length(v1) # Number of elements of v1

cos(v1) #Cosine performed component by component

exp(v1) #Exponential performed component by component

# Nota: Stiamo solo stampando a video i risultati delle operazioni. Non abbiamo creato

nuovi oggetti,

# come potete vedere dall'environment, né modificato v1 o v2

# Per creare o modificare oggetti, bisogna SEMPRE effettuare un "assegnamento",

ovvero usare il

# simbolo =

v1 = 2*v1 # Adesso ho modificato v1

v3 = v1*v2 # Adesso ho creato un nuovo vettore, v3, prodotto di v1 e v2.

#POSITION INDICES

min(data$weight)

max(data$weight)

range(data$weight)

mean(data$weight) # la media

median(data$weight) # mediana

quantile(data$weight, c(0.25, 0.75)) # quantili a 0.25 e 0.75

summary(dati$peso) # Global Visualization

#DISPERSION INDICES

var(data$weight) # sample variance / varianza campionaria

sd(data$weight) #sample standard / deviazione standard campionaria come

IQR(data$weight) # Interquartile Range / Range interquartile

#HISTOGRAMS: it counts the frequency of data in each arbitrary numerical interval.

x11()

hist( data$weight, prob = T ) # on the y-axis there are relative frequencies

abline(v= median(data$weight), col = 'red') # a red vertical line at the median

abline(v= mean(data$weight), col = 'green') # a green vertical line at the mean

#let's specify the number of bins with "breaks = ..."

hist( data$weight, prob = TRUE, breaks = 15, col = 'orange',

main = 'Histogram of weight', xlab = 'Weight', ylab = 'Density' )

#When I don't know what a function does I use help()

help(hist)

#BOXPLOT: useful for identifying Asymmetries or the presence of OUTLIERS

#Upper Limit = Q3 + 1.5*IQR

#Lower Limit = Q3 - 1.5*IQR

x11()

boxplot( data$weight, main = 'Boxplot of weight' , col = "forestgreen")

#Boxplot with the observations of the variable "Haircut"

# for the two subsamples identified by sex.

x11()

boxplot(students$Haircut ~ students$Sex, col=c('pink', 'lightblue'))

#oppure (per averne due vicini)

x11()

par(mfrow = c(1,2)) # sara' un device 1x2, quindi c'e' spazio per 2 grafici.

for( i in 1:2)

{

hist(dati[,i], main = names(dati)[i], xlab = names(dati)[i], prob = T)

# Per riportare le frequenze metto prob = T. Lascio che i bin vengano costruiti

autonomamente

}

#CATEGORICAL VARIABLES

# Absolute frequencies

abs_freq = table(data$sex)

abs_freq

# Relative Frequencies

#method 1

rel_freq = abs_freq / length(data$sex)

rel_freq

#method 2

rel_freq2 = prop.table(abs_freq)

rel_freq2

#Analysis differentiated by class

# Weight of only male penguins:

weight_males = data$weight[data$sex == "M"]

weight_females = data$weight[data$sex == "F"]

# Respective means

mean(weight_males)

mean(weight_females)

# I can also directly create two dataframes and analyze them separately

male_penguins = data[data$sex == "M", ] #we are considering all the columns but

not all the rows, just the ones which respect data$sex == "M

female_penguins = data[data$sex == "F", ]

mean(male_penguins$weight)

#EXERCISE EXAMPLE - Repeatition

prov_data =

factor(c(rep("Milano",60),rep("Lecco",20),rep("Bergamo",18),rep("Brescia",22)))

# Shuffle

prov_data = sample(prov_data)

# These are categorical data

# Identify values

factor(prov_data)

# Calculate absolute frequencies

num = table(prov_data)

# Calculate relative frequencies

table(prov_data)/length(prov_data)

# Bar chart

barplot(table(prov_data)/length(prov_data), col = "gold")

# Pie chart

pie(table(prov_data))

# Calculation of the mode

num[num == max(num)]

#boxplot of the penguins’ weight, divided by species:

boxplot(data$weight~data$species, col = "gold")

#------------------------------------------------------------------------------------------------------------------

#LABORATORIO 2# (esplorazione multivariata, ANOVA)

#GRAPHICAL EXPLORATION OF MULTIVARIATE DATA

x11() #scatterplot base

plot(data) #visualizing the interactions between variables

#the function ggpairs di GGally is more flexible.

#

#è riportato il risultato di un test di correlazione di Pearson fatto fra le variabili relative

ad ogni incrocio. Il test è fatto sia per la popolazione completa (valore in nero), che per

le sottopopolazioni relative a ciascun genere.

#' Il valore numerico rappresenta la stima della correlazione, gli asterischi riportano

invece il livello di significatività del test di correlazione [RECAP stat 1: *** vuol dire

MOLTO significativo].

library(GGally)

library(ggplot2)

x11()

ggpairs(data, aes(alpha = 0.4, color = Sex))

#Visualize the weight data relating them to the Sex and Species factors

x11()

ggpairs(penguins, aes(alpha = 0.4, color = species), columns = 1:2) # look only at

speciess as categorical

ggpairs(penguins, aes(alpha = 0.4, color = sex), columns = c(1,3)) # look only at sex

as categorical

ggpairs(penguins, aes(alpha = 0.4, color = species), columns = 1:3) # visualize both

but color by species

# Calcolo della matrice di covarianza

cov(data[ ,-4]) # RICORDA: con dati[,-4] sto prendendo tutte le righe di dati e tutte le

colonne meno la quarta

# Calcolo della matrice di correlazione

cor(data[,-4])

# NORMALITY TESTING

#In a good number of situations where it is necessary to produce confidence intervals,

or perform hypothesis

#tests, it is required that at least one of the following conditions is satisfied:

# - That the observations are distributed as a normal

# - That the observations are numerous

## Data

chol = data$Cholesterol

mu_hat = mean(chol) #mean

sd_hat = sd(chol) #standard deviation

#First Tool (Graphical) -> histogram of the obsverations and overlaying on it the

normal density with the same mean and variance

t = seq(min(chol), max(chol), length.out = 1000) #observations

f = dnorm(t, mean = mu_hat, sd = sd_hat) #normal density

hist(chol, breaks = 15, col = "gold", probability = T, main = "Cholesterol Histogram")

#histogram

lines(t, f, col = "red", type = "l", lwd = 2) #drawing the gaussian

#-Second Tool (Graphical) -> qqplot with the function qqnorm() for quantiles obtained

from the observations

#function qqline() to plot the line of theorical quantiles of the normal distribution

#overlaying the empirical quantiles obtained from the observations

#with the line of theoretical quantiles of the normal distribution

x11()

qqnorm(chol)

qqline(chol, col = "red") ## In case of perfect normality, the data lie on the line

#-Tool 3 (test) -> Shapiro-Wilk Test -> Quantitative Method

#H0: normality of the data

#H1: no normality

## H0 : X~N vs H1

## Low p-value: H0 rejected -> no normality

## High p-value: NO rejection H0 -> normality

shapiro.test(chol)

## Se non c'è evidenza per rifiutare H0, si può lavorare sotto l'ipotesi di gaussianità

#Valori tipici per α:

#0.05 (più comune)

#0.01 (più severo)

#0.10 (più permissivo)

# HYPOTHESIS TESTING ON THE MEAN OF A SINGLE

POPULATION-----------------------------------------------------------

#The test statistic follows a t distribution with n − 1 degrees of freedom.

# - p-value ALTO: l'evidenza empirica non è sufficientemente contraria all'ipotesi nulla

H0, che quindi non può essere rifiutata

# - p-value BASSO: l'evidenza empirica è fortemente contraria all'ipotesi H0 che quindi

va rifiutata

# CASO 1: One Small Sample from a Normal Population — t-test (BEFORE DOING THE

TEST THE NORMALITY ASSUMPTION MUST BE VERIFIED)

#MANUALLY

#data

mu_0 <- 100

n <- length(f_glucose)

mu_hat <- mean(f_glucose)

sd_hat <- sd(f_glucose)

t_stat <- (mu_hat - mu_0) / (sd_hat / sqrt(n))

t_stat

#Left-Tailed Test: �0 ∶ � ≥ 100 vs �1 ∶ � < 100

p_value <- pt(t_stat, df = n - 1)

p_value

#Right-Tailed Test: �0 ∶ � ≤ 100 vs �1 ∶ � > 100:

pval <- 1 - pt(t_stat, df = n - 1)

pval

#Two_Tailed test: �0 ∶ � = 100 vs �1 ∶ � ≠ 100:

pval <- 2 * (1 - pt(abs(t_stat), df = n - 1))

pval

#Two_Tailed Test: �0 ∶ � = 100 vs �1 ∶ � ≠ 100:

#Function

t.test(dati$Colesterolo, mu = 200, alternative = "less") # H0: MU >= X UNILATERO A

SINISTRA - LEFT-TAILED

t.test(dati$Colesterolo, mu = 200, alternative = "greater") # H0: MU <= X UNILATERO

A DESTRA -RIGHT-TAILED

t.test(dati$Colesterolo, mu = 200, alternative = "two.sided") # H0 != 0 BILATERO

-BILATERAL

# - CONFIDENCE INTERVAL FOR THE MEAN

alpha = 0.1 #90% Confidence Interval

quantile = qt(1-alpha/2, df = n-1) # qui utilizzo i quantili (qt anzichè pt)

c(mu_hat - quantile / sqrt(n) * sd_hat,

mu_hat + quantile / sqrt(n) * sd_hat)

# CASO 2: Large Sample from a Normal Population — z-test

#normality assumption verified

#population is large we can do the Z-test

#LEFT-TAILED TEST

#�0 ∶ � ≥ 200 vs �1 ∶ � < 200

mu_0 <- 200

n <- length(data$Cholesterol)

mu_hat <- mean(data$Cholesterol)

sd_hat <- sd(data$Cholesterol)

t_stat <- (mu_hat - mu_0) / sqrt(sd_hat^2 / n)

pval <- pnorm(t_stat)

pval

#RIGHT-TAILED TEST

#�0 ∶ � ≤ 200 vs �1 ∶ � > 200:

mu_0 <- 200

n <- length(data$Cholesterol)

mu_hat <- mean(data$Cholesterol)

sd_hat <- sd(data$Cholesterol)

t_stat <- (mu_hat - mu_0) / sqrt(sd_hat^2 / n)

pval <- 1 - pnorm(t_stat)

pval

#P-VALUE TEST

pval <- pt(t_stat, df = 48) #LEFT TAILED

pval <- 1 - pt(t_stat, df = 48) #RIGHT TAILED

pval <- 2*(1 - pt(t_stat, df = 48)) #TWO SIDED

pval

## BILATERAL TEST

## H0: mu=157 vs H1: mu!=157

mu_0 = 157

n = length(dati$Peso)

mu_hat = mean(dati$Peso)

sd_hat = sd(dati$Peso)

t_stat <- (mu_hat - mu_0) / sqrt(sd_hat^2 / n) #Statistica test

pval = 2*(1-pnorm(abs(t))) #Calcolo il pvalue del test bilatero

pval

#CONFIDENCE INTERVAL FOR THE MEAN

alpha <- 0.1

quantile <- qnorm(1 - alpha / 2)

c(mu_hat - quantile / sqrt(n) * sd_hat,

mu_hat + quantile / sqrt(n) * sd_hat)

#CASO 3: NON - NORMAL POPULATION WITH LARGE SAMPLE - ASYMPTOTIC Z-TEST

#Not Normally Distributed but since the sample size is large, the Central Limit

Theorem guarantees

#that the sample mean is approximately Gaussian — and we can therefore proceed

with an asymptotic

#z-test.

#TWO-TAILED TEST

#�0 ∶ � = 157 vs �1 ∶ � ≠ 157

mu_0 <- 157

n <- length(data$Weight)

mu_hat <- mean(data$Weight)

sd_hat <- sd(data$Weight)

t_stat <- (mu_hat - mu_0) * sqrt(n) / sd_hat

pval <- 2 * (1 - pnorm(abs(t_stat)))

pval

#Right-tailed and left-tailed variants follow the same logic seen in the previous

examples.

##### CONFIDENCE INTERVAL FOR THE MEAN (95%)

alpha <- 0.05

quantile <- qnorm(1 - alpha / 2)

n <- length(data$Weight)

c(mu_hat - quantile / sqrt(n) * sd_hat,

mu_hat + quantile / sqrt(n) * sd_hat)

# 2 POPULATION - COMPARISON TEST BETWEEN MEANS OF INDEPENDENT NORMAL

POPULATIONS

#�0 ∶ �� = �� vs �1 ∶ �� ≠ ��

#dataset divided into 2 groups

cholesterol_men <- data$Cholesterol[data$Sex == "M"]

cholesterol_women <- data$Cholesterol[data$Sex == "F"]

#ASSUMPTION 1: NORMALITY

shapiro.test(cholesterol_men)

shapiro.test(cholesterol_women)

# - CASO 1 : Unkown Variances but equal (to be verified)-> sample variances

#ASSUMPTION 1 -> FISHER'S TEST for homogeneity of variances (var.test())

#H0: homogeneity of variances

#H1: different variances

#low p-value -> rejection -> no homogeneity of variances

var.test(cholesterol_men, cholesterol_women)

#formula

t.test(cholesterol_men, cholesterol_women, var.equal = TRUE)

t.test(cholesterol_men, cholesterol_women, var.equal = TRUE, alternative = "less")

t.test(cholesterol_men, cholesterol_women, var.equal = TRUE, alternative = "greater")

t.test(cholesterol_men, cholesterol_women, var.equal = TRUE, alternative =

"two.sided")

#manually

cholesterol_men <- data$Cholesterol[data$Sex == "M"]

cholesterol_women <- data$Cholesterol[data$Sex == "F"]

mean_diff <- mean(cholesterol_men) - mean(cholesterol_women)

n_men <- length(cholesterol_men)

n_women <- length(cholesterol_women)

var_men <- var(cholesterol_men)

var_women <- var(cholesterol_women)

var_pooled <- ((n_men - 1) * var_men +

(n_women - 1) * var_women) / (n_men + n_women - 2)

sd_pooled <- sqrt(var_pooled)

t <- mean_diff / (sd_pooled * sqrt(1/n_men + 1/n_women))

pvalue <- 2 * (1 - pt(abs(t), df = n_men + n_women - 2))

pvalue

df = m_men + n_woman - 2

t_critico = qt(0.99, df = df)

t_critico

#NOTA:

#1) If variances are unknown and unequal we do the test with var.equal = FALSE

t.test(cholesterol_men, cholesterol_women, var.equal = FALSE)

#2) If population is not normally distributed we conduct a Z-test instaed of T-test

#example z-test

#H0 = weight male <= weight female

#' H1 = male weight > female weight

#' # Let's see it graphically again

ggpairs(penguins, aes(alpha = 0.4, color = sex), columns = c(1, 3)) # look only at sex

as categorical

tapply(penguins$weight, penguins$sex, mean)

# visually the means seem very similar

men <- penguins$weight[penguins$sex == "M"]

women <- penguins$weight[penguins$sex == "F"]

#z-test

# Manually construct the z-test

X = females$weight

Y = males$weight

nX = length(X)

nY = length(Y)

s2X = var(X)

s2Y = var(Y)

s2p = 1 / (nX + nY - 2) * ((nX - 1) * s2X + (nY - 1) * s2Y)

sp = sqrt(s2p)

z = (mean(Y) - mean(X)) / sqrt(s2X / nX + s2Y / nY)

z

pval = 1 - pnorm(z)

pval

# - CASO 2 : 2 POPULATION - Known variances

#no formula, only manually

#both population are normally distributed

mean_diff <- mean(cholesterol_men) - mean(cholesterol_women)

n_men <- length(cholesterol_men)

n_women <- length(cholesterol_women)

# Suppose the variances are KNOWN and we denote them by var_M and var_F

var_M <- 150

var_F <- 110

z <- mean_diff / sqrt(var_M / n_men + var_F / n_women)

pvalue <- 2 * (1 - pnorm(abs(z)))

pvalue

#NOTA: In this case as well, one-sided tests and confidence intervals for the difference

of means can be const