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)

## alt.: leggi i dati da file di testo (qui: comma separated values)
Tab6.1 <- read.csv("Tab6.1.csv")

## 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)
## vedi anche esercizio Cap6.R

## Abbiamo stabilito che la correlazione è significativa.
## Ipotizziamo una relazione lineare e stimiamo il modello
## Costo = alfa  + beta*Produzione + u

## la produzione è espressa in centinaia (per avere numeri
## "più belli da leggere"); la funzione I() "protegge" il
## calcolo da interpretazioni sbagliate.

## specifico il modello creando una formula
fm <- Costo ~ I(Produzione/100)
## stimo il modello a OLS con la funzione lm()
mod <- lm(fm, Tab6.1)
## ispeziono i risultati
summary(mod)

## coefficienti:
coef(mod)

## intervallo di confidenza per beta:

## valore puntuale:
beta.hat <- coef(mod)[2]
## ES(beta):
SE.beta <- summary(mod)$coef[2,2]

## (valore assoluto del-) valore critico al 95%, due code:
z.bar <- qt(0.975, df=n-2)
## estremi del C.I.
ci.lower <- beta.hat - z.bar*SE.beta
ci.higher <- beta.hat + z.bar*SE.beta

## rappresentazione grafica:
plot(Tab6.1$Produzione/100, Tab6.1$Costo, pch=19, col="red",
     xlim=c(0, 46), ylim=c(0, 36))
abline(mod, lwd=2)

## Analisi grafica dei residui:
qqnorm(resid(mod)) # nota: il grafico di R inverte gli assi
plot(resid(mod), col="blue", type="l")
abline(h=0, lty=2, col="orange")