Copiare i comandi che seguono nella Console di R
(versione del 12/03/2003)
- Distribuzioni di frequenza
load("dati1.rda")
str(dati)
attach(dati)
table(X) # frequenze assolute
table(X)/length(X) # frequenze relative
table(X)/length(X)*100 # frequenze percentuali
cumsum(table(X)) # frequenze cumulate
# analogamente per Y e Z
table(W) # non molto efficace
table( cut(W, breaks=c(40,50,58,70,95)) ) # dati raccolti in classi
classi <- c(40,50,58,70,95) # definiamo i punti estremi delle classi
table( cut(W, br=classi, right=FALSE) ) # classi chiuse a sinistra
hist(W, br=classi, plot = FALSE) # informazioni aggiuntive comprese le densita' di frequenza
table(X,Y) # distribuzioni doppie
table(X, cut(W,br=classi))
table(X,Y,Z) # tabelle a piu' vie
ftable(X,Y,Z) # flat tables
tab <- table(X,Y)
tab # distribuzione congiunta X su Y
margin.table(tab,1) # marginale di X
margin.table(tab,2) # marginale di Y
prop.table(tab) # congiunte relative
prop.table(tab,1) # condizionate di Y a X (relative)
prop.table(tab,2) # condizionate di X a Y (relative)
Rappresentazioni grafiche
# in alternativa un metodo piu' efficiente e' "hist"
hist(W, br=classi) # disegno dell'istogramma
hist(W, br=classi, freq=TRUE) # forziamo il disegno dell'istogramma sbagliato! R ci mettera' in guardia
hist( W, br=classi ) # classi chiuse a sinistra
hist( W, br=3 ) # tre classi di stessa ampiezza
hist( W, br="ST" ) # classi ottenute con il metodo di Sturges
plot(X) # diagramma a barre
plot(table(Z)) # diagramma a bastoncini
pie(table(Y)) # diagrammi a torta
pie(table(Y), density=15, angle=15+10*(1:4), col=1:4) # diagrammi a torta (2)
plot(Z) # grafico assurdo, Z e' quantitativo
plot(Z,W) # diagrammi di dispersione
plot(X,Y) # aste impilate
plot(table(X,Y)) # rappresentazioni di tabelle di contingenza
boxplot(W) # boxplot
Poligono di frequenza
hist.pf <- function(x, br, ...){ # i puntini "..." indicano parametri generici che vengono passati ad altre funzioni
if(missing(br)) # manca il valore di "br"
ist <- hist(x,...)
else
ist <- hist(x, breaks=br,...)
if(ist$equidist)
lines( c(min(ist$breaks),ist$mids,max(ist$breaks)),
c(0,ist$counts,0))
else
lines( c(min(ist$breaks),ist$mids,max(ist$breaks)),
c(0,ist$intensities,0))
}
hist.pf(W)
library(labstatR) # R ci segnala che hist.pf esiste anche nella libreria labstatR
plot(density(W)) # stima kernel della densita' (bandwidth ottimale)
plot(density(W,bw=3)) # stima kernel della densita' (bandwidth imposto)
density(W)
hist.pf(W)
lines(density(W),col="red")
bubbleplot(table(X,Y)) # bubbleplot a la Excel
Funzione di ripartizione empirica
library(stepfun)
plot(ecdf(Z)) # ecdf: Empirical Cumulative Distribution Function
ecdf(Z)(2) # e' una vera e propria funzione quindi possiamo calcolarla nel punto 2
# disegnamola per dati continui
classi <- c(30, 40, 50, 58, 70, 95, 100)
Fi <- cumsum( table( cut(W,classi) ) ) / length(W)
Fi <- c(0, Fi)
plot(classi, Fi, type = "b", axes = FALSE, main = "Funzione di ripartizione")
axis(2, Fi)
axis(1, classi)
box()
Q-Q plot
qqnorm(W) # quantili di W contro quelli della Gaussiana
qqline(W)
qqplot(W,Z) # Q-Q plot W contro Z
Statistiche di base
summary(X) # distribuzione di frequenza
summary(Z) # quantili, media, range
summary(table(X,Y)) # test del Chi-Quadrato di indipendenza
var(W) # varianza (corretta) di W
cor(Z,W) # coefficiente di correlazione
cbind(Z,W,X) # crea una matrice accostando le colonne, trasforma X in numerico
cor( cbind(Z,W,X) ) # matrice di correlazione (Attenzione: qui X e' trattato come numerico!)
Y
sort(Y) # Y e' ordinale ed R sa come trattarlo
Me(Y) # calcola la mediana per fenomeni qualitativi (labstatR)
median(W) # calcolo della mediana
quantile(W, c(0.1,0.3,0.7,0.93,.98)) # calcolo dei quantili
detach(dati)
rm(list=ls())
Indice di concentrazione e spezzata di Lorenz
x <- c(1,1,1,4,4,5,7,10)
y <- c(1,1,1,1,1,4,4,4,5,9,100,100,200)
gini(x,col="blue") # calcolo dell'indice di concentrazione di Gini e disegno della curva di Lorenz
gini(y,add=TRUE,col="red") # stessa cosa per "y" sovrapponendo alla precedente la nuova curva di concentrazione
Verifica di ipotesi e intervalli di confidenza
t.test(1:10, y = 7:20) # t-test per due campioni, liv confidenza 95%, alternativa medie diverse
t.test(1:10, y = 7:20, alt="gr") # come sopra ma con alternativa > 0
t.test(1:10, y = 7:20, conf=0.99) # liv confidenza 0.99
heads <- rbinom(1, size = 100, pr = 0.5)
prop.test(heads, 100) # basato sul test ChiQuadrato
Modello lineare
data(cars)
str(cars)
plot(cars)
lm(dist ~ speed, data=cars) # Linear Model (lm) dist = f(speed), "dist" e "speed" nel dataset "cars"
model <- lm(dist ~ speed, data=cars)
summary(model)
plot(model)
str(model) # struttura complessa di un modello lineare
plot(model$res) # grafico dei residui
plot(cars)
abline(model) # grafico della retta di regressione
predict(model) # valori stimati dalla retta di regressione
predict(model, newdata=data.frame(speed=9)) # valore del modello per speed=9
predict(model, newdata=data.frame(speed=c(9,2,13))) # stima in corrispondenza dei valori di speed 9,2,13
pc <- predict(model,interval="c") # calcolo delle bande di confidenza
pc
attach(cars)
plot(speed,dist)
abline(model)
matlines(speed,pc[,2:3],lty=3:3,col="red",lwd=2:2) # disegno delle bande di confidenza
detach(cars)
library(MASS)
data(whiteside)
help(whiteside)
coplot(Gas ~ Temp | Insul,whiteside) # regressione di Gas su Temp condizionatamente a Insul
gasB <- lm(Gas ~ Temp, whiteside, subset = Insul=="Before")
gasA <- lm(Gas ~ Temp, whiteside, subset = Insul=="After")
summary(gasA)
coef(gasA) # fornisce i coefficienti
resid(gasA) # estrae i residui dal modello
hist(resid(gasA)) # istogramma dei residui
fitted(gasA) # estrae i valori stimati
gA <- subset(whiteside,Insul=="After")
gB <- subset(whiteside,Insul=="Before")
pAdatoB <- predict(gasB, newdata=gA) # predice con il modello "Before" i valori per il sottogruppo "After"
plot(gA[,2:3])
points(gB[,2:3], pch=18, col="red")
Cambiamenti di scala
attach(cars)
plot(speed,dist)
cor(speed,dist)
plot(speed,sqrt(dist))
cor(speed,sqrt(dist))
summary(lm(dist ~ speed))
summary(lm(sqrt(dist) ~ speed))
Regressione non parametrica
library(modreg) # modern regression
plot(speed, dist)
lines(ksmooth(speed, dist, "normal", bandwidth=2)) # Nadaraya-Watson (kernel)
lines(ksmooth(speed, dist, "normal", bandwidth=5),lty=3)
lines(ksmooth(speed, dist, "normal", bandwidth=10))
plot(cars)
lines(lowess(cars)) # regressione locale robusta
lines(lowess(cars, f=.2), lty = 3)
legend(5, 120, c(paste("f = ", c("2/3", ".2"))), lty = c(1:3))
Anova a 1 via
x <- c(rep(1,6), rep(2,6), rep(3,6), rep(4,6))
y <- c(64, 72, 68, 77, 56, 95, 78, 91, 97, 82, 85, 77,
75, 93, 78, 71, 63, 76, 55, 66, 49, 64, 70, 68)
bombe <- data.frame(y,x)
lm(y~x)
anova(lm(y~x))
plot(x,y)
abline(lm(y~x))
x <- factor(x)
x
lm(y~x)
anova(lm(y~x))
Specificazioni dei modelli
x1 <- runif(50) # generatore di numeri casuali
x2 <- sin(runif(50))
x3 <- rnorm(50,sd=45)
x4 <- rchisq(50,df=1)
x5 <- rt(50,df=3)
y <- 4 + 2*x1 - x2 + rnorm(50)
lm(y ~ x1 + x2)
z <- data.frame(y, x1, x2, x3, x4, x5)
lm( y ~ . , data = z)
# specificando i regressori da includere
lm(y ~ x1 + x3 + x5, data=z)
# o specificando quelli da escludere
lm(y ~ . - x2 - x4, data=z)
lm(y ~ x1 + x2 -1, data=z)
Modelli lineari generalizzati: regressione logistica
y <- y < mean(y) # trasformiamo y in una successione di 0 e 1
glm(y~x1+x2+x3, family=binomial("logit"))
glm(y~x1+x2+x3, binomial)
data(interinale)
str(interinale,vec.len=0)
glm(avviato~., binomial, data=interinale) -> model
summary(model) # livello di istruzione e' non significativa
glm(avviato~. - istruzione, binomial, data=interinale) -> model1
summary(model1) # il modello migliora leggermente (AIC piu' basso)
predict(model1,type="response") -> pr
plot(density(pr),xlim=c(0,0.2),main="",xlab="probabilita' di avviamento")
Generazioni di numeri casuali
runif(10) # Uniforme
1*(runif(10) < 0.3) # Ber(0.3)
rbinom(10, 20, 0.4) # Bin(20,0.4)
rnorm(10, mean=5, sd=3) # N(5,3)
rchisq(10, df=2) # Chi2(2)
rt(10, df=3) # t(3)
Tavole statistiche e quantili
pnorm(.3, mean=.1, sd=2) # P(X<.3) X ~ N(.1,4)
dnorm(.4) # densita' di N(0,1) calcolata in .4
qnorm(0.96) # 96-esimo percentile di una N(0,1)
Campionamento e generatori ad hoc di v.c. discrete
x <- 1:10
sample(x,3) # campionamento senza ripetizione
sample(x,11) # errore!
sample(x,11,rep=TRUE) # campionamento con ripetizione
x <- 1:3
sample(x,11,rep=TRUE,prob=c(0.1,0.6,0.3)) # variabile casuale X t.c. P(X=1)=0.1, P(X=2)=0.6, P(X=3)=0.3