Statistica III - Pratica (comandi di R)

Come importare file su R

Impostazione della directory di lavoro

Session > set working directory > choose directory. NB: Devo aver salvato il file come .csv, .txt, .xls, .xlsx nella directory che vado a scegliere.

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 <- read.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 test
y.i <- Y/n.i #frequenze relative dei successi

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

Diagramma di dispersione dei valori delle freq. relative osservate vs temperatura

plot(y.i~X + X)

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

ng <- 3 #numero dei gruppi che mi fornisce il test
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])

Funzione che per ogni gruppo fornisce S:

obs.exp <- function(g, num.g) {
  succ.obs <- sum(g[,1])
  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, breaks=6)

Dichotomizzare 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

Individuare 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

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

Fornire tabella a doppia entrata che riporta Y distinta per X e X

tab <- xtabs(Y ~ X + X, data=data)
prop.table(tab) #tabella a doppia entrata frequenze relative

Boxplot che riporti andamento della Y condizionato per X e X

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

Boxplot per X distinta per le due modalità assumibili dalla X

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

Stimare un opportuno GLM per Poisson

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

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₂.

m <- nrow(data)
n <- sum(Y1 + Y2)
n.i <- Y1 + Y2
y.i <- Y1 / n1
mod <- glm(cbind(Y1, Y2) ~ X + X + X, family=binomial, data=data)
mod <- glm(y.i ~ X + X + X, family=binomial, weights=n.i, data=data) #metodo migliore

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 bontà 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 corrispondente.

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 è 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 è distribuita come una Poisson

Distribuzione empirica:

tab <- xtabs(~Y, data=data)
ascisse <- as.numeric(names(tab))
plot(ascisse, tab)
points(ascisse, tab) #dist. empirica