library(XLConnect)
wb <- loadWorkbook("~/Dropbox/StatAz/Esercizi/TabelleVuote.xlsx") 

## Leggi la tabella dei dati (crea un oggetto data.frame)
Tab6.1 <- readWorksheet(wb, sheet = "Cap6", region="A3:C25",
                    header = TRUE)


## Ispezione grafica della correlazione tra produzione e costo
plot(Tab6.1$Produzione, Tab6.1$Costo, pch=19, col="green3")
## (su una scala "neutrale":)
plot(Tab6.1$Produzione, Tab6.1$Costo, pch=19, col="red",
     xlim=c(0, 4600), ylim=c(0, 36))

## calcolo del coefficente di correlazione (campionario):
n <- dim(Tab6.1)[[1]]

## rinominiamo solo per comodità:
x <- Tab6.1$Produzione
y <- Tab6.1$Costo

## covarianza (campionaria):
cov.xy <- sum((x-mean(x))*(y-mean(y)))/(n-1)
## automaticamente, in R:
cov(x,y)

## varianze (campionarie):
var.x <- sum((x-mean(x))^2)/(n-1)
var.y <- sum((y-mean(y))^2)/(n-1)

## coefficiente di correlazione (campionario, ma la formula è la stessa)
r.xy <- cov.xy/(sqrt(var.x)*sqrt(var.y))

## Inferenza sul coefficiente di correlazione: 
## stimiamo il suo errore standard sotto l'ipotesi H_0: r=0
ES.r <- sqrt((1-r.xy^2)/(n-2))

## da cui la statistica t=r.xy/ES.r per la medesima ipotesi:
t.r <- r.xy*sqrt((n-2)/(1-r.xy^2))

## rifiutare o non rifiutare?
## calcolo i valori critici per n-2 gradi di libertà:
tcrit <- qt(0.975, df=n-2)
## ... e confronto:
abs(t.r) > abs(tcrit)

## graficamente:
plot(density(rt(1000000, df=n-2)), col="blue", xlim=c(-10, 10),
     main="Distribuzione di t.xy sub H0: rho=0", xlab="", ylab="")
## valori critici:
abline(v=-tcrit, col="red")
abline(v=tcrit, col="red")
## dove cade la statistica test
abline(v=t.r, col="green3")


Tab6.2 <- readWorksheet(wb, sheet = "Cap6", region="A29:L51",
                        header = TRUE)

## Ispezione grafica simultanea della correlazione tra numerose variabili
pairs(Tab6.2[,-1]) # elimino la colonna dei nomi (la prima)
library(corrplot)
corrplot(cor(Tab6.2[,-1]))



Tab6.6 <- readWorksheet(wb, sheet = "Cap6", region="A55:D89",
                        header = TRUE)