Statistica III - Pratica (comandi di R)

COME IMPORTARE FILE SU R

1. Session > set working directory > choose directory

NB: devo aver salvato il file come .csv .txt .xls . xlsx nella directory che vado a scegliere

2. Impostare il percorso e la directory direttamente con i comandi R:

setwd("percorso completo del file")

data <- read.csv (“file.csv”, header =TRUE, sep=”;”)

data <- raed.table(“file.txt”, header=TRUE, sep=”;”)

data <- read.table(“file.xls”, header=TRUE, sep=”;”)

str(data)

attach(data)

View(data)

GLM PER DATI SPARSI

m <- nrow(data) #numero valori esplicativa

n.i <- rep(6,m) #6=n° prove per ogni X, fornito dal testo

y.i <- Y/n.i #frequenze relative dei successi

mod <- glm(y.i~ X +X , family=binomial, weights=n.i, data=data)

1 2

summary(mod)

#diagramma di dispersione dei valori delle freq. relative osservate vs temperatura

plot(y.i~X +X )

1 2

RIPORTARE FREQUENZE ASSOLUTE OSSERVATE E ATTESE/PREVISTE

attese <- n.i*predict(mod, type="response")

frequenze <- cbind(Y, attese.logit,)

TEST HOSMER-LEMESHOW PER 3 GRUPPI 8,8,7=numerosità (da usare in caso di dati sparsi)

H : il modello è buon interprete dei dati

0

ng <- 3 #numero dei gruppi che mi fornisce il testo

g1 <- frequenze[1:8,] ; num.g1<-sum(n.i[1:8])

g2 <- frequenze[9:16,] ; num.g2<-sum(n.i[9:16])

g3 <- frequenze[17:23,] ; num.g3<-sum(n.i[17:23])

obs.exp <- function(g,num.g){ #funzione che per ogni gruppo fornisce S ,

obs

succ.obs <- sum(g[,1]) S , I , I per poter calcolare X.HL

exp obs exp

succ.exp <- sum(g[,2])

ins.obs <- num.g-succ.obs

ins.exp <- num.g-succ.exp

c(succ.obs,succ.exp,ins.obs,ins.exp)

}

#matrice che per ogni gruppo contiene i 4 valori

matrice <-cbind(obs.exp(g1,num.g1),obs.exp(g2,num.g2),obs.exp(g3,num.g3))

X.HL <- sum((matrice[1,]-matrice[2,])^2/matrice[2,])+sum((matrice[3,]matrice[4,])^2/matrice[4,])

pvalue <- 1-pchisq(X.HL,g-2)

VALUTARE SE CI SONO VALORI NA E CONSIDERO SOLO OSSERVAZIONI COMPLETE

anyNA(data)

dataC<- data [complete.cases(data),]

anyNA(dataC)

RICLASSIFICARE UNA VARIABILE IN 6 CLASSI

data$variabile<- cut(variabile, beaks=6)

DICOTOMIZZARE UNA VARIABILE

data$var <- factor (ifelse (data$var < soglia , 0 , 1), labels=c(“scarsa” , ”elevata”))

#in questo modo il dataset aggiunge una colonna chiamata var

data<-data[,c(1,2,5,7)] #in questo modo ora posso lavorare con un nuovo dataset che

ha solo le variabili relative alle colonne che ho selezionato

INDIVIDUO POSIZIONE DELLE OSS. CHE NON HANNO MAI VIAGGIATO(SERVICE) E

LE ELIMINO DAL DATASET

pos.0 <- which(data$service==0)

data <- data[-pos.0,] #elimino le navi che non hanno mai viaggiato

summary(data)

COSTRUIRE UN INTERACTION PLOT CHE RIPORTI ANDAMENTO DI Y AL VARIARE DI

X DISTINTO PER IL VALORE DI X

1 2

with(data, interaction.plot(X ,X , Y))

1 2

FORNIRE TABELLA A DOPPIA ENTRATA CHE RIPORTA Y DISTINTA PER X E X

1 2

tab <- xtabs(Y~X +X , data=data)

1 2

prop.table(tab) #tabella a doppia entrata frequenze relative

BOXPLOT CHE RIPORTI ANDAMENTO DELLA Y CONDIZIONATO PER X E X

1 2

boxplot(Y~X +X , data=data)

1 2

BOXPLOT PER X DISTINTA PER LE DUE MODALITÀ ASSUMIBILI DALLA X

1 2

plot(factor(X ), X ) #se la variabile é factor, il plot riporta il boxplot

2 1

STIMARE UN OPPORTUNO GLM PER POISSON

mod <- glm(Y~offset(log(X )) +X +X , family=poisson, data=data)

1 2 3

#uso il + se SENZA INTERAZIONE, uso il * se HO INTERAZIONE

#offset SOLO SE c’è una variabile che rappresenta il n° delle esposizioni

à

STIMARE UN OPPORTUNO GLM PER BINOMIALE Y è divisa in Y e Y

1 2

m <- nrow(data)

n <- sum(Y +Y )

1 2

n.i <- Y +Y

1 2

y.i <- Y /n

1

mod <- glm (cbind(Y ,Y ) ~ X +X +X , family=binomial, data=data)

1 2 1 2 3

mod <- glm(y.i ~ X +X +X , family=binomial, weights=n.i, data=data) #metodo migliore

1 2 3

#uso il + se SENZA INTERAZIONE, uso il * se HO INTERAZIONE

#se non specifico il link in automatico viene considerato quello canonico

mod2 <- update(mod, subset = (Dept != "A")) #aggiorno mod togliendo oss. dove Dept=A

TEST PER LA BONTA’ DEL MODELLO BASATO SULLA DEVIANZA

dev <- mod$deviance

gdl <- mod$df.residual

p.value <- 1-pchisq(dev, gdl)

QUALE PATTERN CONTRIBUISCE MAGGIORMENTE ALLA DEVIANZA DEL

MODELLO? Riporto sia il pattern che il valore d corrispondente

i

pattern.max <- which.max(residuals(mod)^2)

residuals(mod)[pattern.max]^2

(pattern.max/dev) #% della devianza che rappresenta

DECIDERE TRA DUE MODELLI CON 2 TEST DIVERSI

anova(mod, mod2, test='Chisq') #confronto con il modello precedente

à

#il p-value è quasi zero rifiutiamo il modello

precedente in favore di quello nuovo

f.obs <- (dev-dev2)/(dev2/gdl2) #alternativa all’ANOVA ma è non precisa

1-pf(f.obs,1,gdl2)

AIC(mod2) #per il criterio AIC scegliamo il modello che ha AIC minore

AIC(mod3)

DECIDERE SE UNA VARIABILE RISPOSTA E’ NORMALMENTE DISTRIBUITA

shapiro.test (Y) #test

qqnorm (Y) #graficamente

curve(dnorm(x, mean(Y), sd(Y), add=T)

plot (density(Y)) #motivazione teorica

DECIDERE SE UNA VARIABILE RISPOSTA E’ DISTRIBUITA COME UNA POISSON

#distribuzione empirica

tab <- xtabs(~Y, data=data)

ascisse <- as.numeric(names(tab))

plot(ascisse, tab)

points(ascisse, tab) #dist. empirica = di