# Analisi Bivariata - Indici ----

## Due variabili categoriali ----
library(insuranceData)
data(AutoBi)

# sistemiamo le variabili per l'analisi
AutoBi$LOSSclass <- cut(AutoBi$LOSS, breaks = c(0, 0.5, 2, 4, 8, 1100))
AutoBi$ATTORNEY  <- factor(AutoBi$ATTORNEY)
levels(AutoBi$ATTORNEY) <- c("yes", "no")
AutoBi$CLMSEX    <- factor(AutoBi$CLMSEX)
levels(AutoBi$CLMSEX)   <- c("M", "F")

# tabella con la distribuzione congiunta di ATTORNEY e CLMSEX
tab1 <- table(AutoBi$ATTORNEY, AutoBi$CLMSEX)
tab1

# distribuzioni condizionate per riga e colonna
rtab <- prop.table(tab1, margin = 1)  
ctab <- prop.table(tab1, margin = 2)  
rtab; ctab

# visualizziamo le distribuzioni condizionate per colonna
par(mfrow = c(1, 2))
barplot(ctab, legend = TRUE, beside = TRUE, ylim = c(0, 1),
        main = "Avvocato | Sesso")

barplot(prop.table(table(AutoBi$ATTORNEY, AutoBi$LOSSclass), margin = 2),
        beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(horiz = T, x = "top"),
        main = "Avvocato | Classe di danno")

par(mfrow = c(1, 1))

# calcolo indice chi quadro X2 a mano
tot  <- sum(tab1)
mr   <- margin.table(tab1, 1)        # marginale di riga
mc   <- margin.table(tab1, 2)        # marginale di colonna
hatn <- mr %*% t(mc) / tot           # frequenze teoriche di indipendenza 
hatn

X2 <- sum((tab1 - hatn)^2 / hatn)
X2

# indice X2 normalizzato
X2.norm <- X2/(tot*min(length(mr)-1, length(mc)-1))
X2.norm

# Equivalente con chisq.test()
chisq.test(tab1, correct = FALSE)

# Test per avvocato vs classe di danno
tab2 <- table(AutoBi$ATTORNEY, AutoBi$LOSSclass)
chisq.test(tab2, correct = FALSE)

X2.LOSS <- chisq.test(tab2, correct = FALSE)$statistic
X2.norm <- X2.LOSS/(tot*min(nrow(tab2)-1, ncol(tab2)-1))
X2.norm


### Esercizio 1 - Titanic ----
data("Titanic")
tab <- margin.table(Titanic, margin = c(1, 4))


## Una variabile continua e una categoriale ----

# visualizzazione tramite boxplot condizionati
data(iris)
boxplot(Sepal.Length ~ Species, data = iris,
        col = c("lightblue", "blue", "darkblue"))
points(1:3, tapply(iris$Sepal.Length,iris$Species, mean), pch = 7, col = "red")

# test ANOVA
testanova <- aov(Sepal.Length ~ Species, data = iris)
summary(testanova)

# calcolo indice eta2
n      <- nrow(iris)
Devtot <- var(iris$Sepal.Length) * (n - 1)
Devtot   # deve coincidere con la somma delle due Sum Sq

DevBeetw <- summary(testanova)[[1]]["Species", "Sum Sq"] # Devianza between

eta2 <- DevBeetw/Devtot  # eta quadro
eta2 

## Due variabili continue ----
library(MASS)
data(whiteside)

# scatterplot
plot(whiteside$Temp, whiteside$Gas, pch = 19,
     xlab = "Temperatura (°C)", ylab = "Consumo di gas (mc)",
     main = "Consumo di gas vs temperatura esterna")


# covarianza
n <- length(whiteside$Temp)

# covarianza descrittiva 
codev  <- sum((whiteside$Temp - mean(whiteside$Temp)) * 
                (whiteside$Gas - mean(whiteside$Gas)))
covar  <- codev / n
covar

# covarianza in R
cov(whiteside$Temp, whiteside$Gas)
codev / (n - 1) 

# Correlazione
cor(whiteside$Temp, whiteside$Gas)