Statistica III - Teoria + Pratica

Come importare file Excel su R

Apri il file Excel, abilita la modifica, salva con formato CSV delimitato da virgole, apri il nuovo file con blocco note, copialo e incollalo in un nuovo file txt e apri quello con R. Oppure imposti la directory da R, e la copi al posto di percorso:

data=read.csv(file='percorso/data.csv', header=T, sep=',')

Decidere se una variabile risposta è normalmente distribuita

Per effettuare il test di normalità utilizza i seguenti comandi:

shapiro.test(data$variabile) # test
qqnorm(data$variabile) # graficamente
plot(density(data$variabile)) # altro strumento grafico

Decidere se una variabile risposta è distribuita come una Poisson

Per valutare se la variabile segue una distribuzione di Poisson, usa questi comandi:

tab=xtabs(~y)
ks.test(jitter(tab), "ppois", mean(y)) # test
qqPlot(tab, distribution="pois", lambda=mean(y)) # graficamente. I punti fuori non sono compatibili con una distribuzione Poisson
par(mfrow=c(1,2)) # altro strumento grafico
ascisse=as.numeric(names(tab))
plot(ascisse, tab, type="h") # distribuzione empirica
teor=rpois(250, mean(y))
tab.teor=xtabs(~teor)
ascisse.teor=as.numeric(names(tab.teor))
plot(ascisse.teor, tab.teor, type="h") # distribuzione teorica

Ricodificare una variabile come qualitativa

data$variabile=as.factor(data$variabile)
levels(data$variabile)=c('nome1','nome2')

Dichotomizzare una variabile

classi=as.factor(ifelse(data$variabile<numero,0,1))
levels(classi)=c('nomeclasse0','nomeclasse1')

Tabella a doppia entrata frequenze assolute

Creare una tabella di contingenza:

tab=xtabs(~x+y)

Tabella a doppia entrata frequenze relative

Creare la tabella delle frequenze relative:

prop.table(tab)
prop.table(tab, 1) # condizionate (1: righe, 2: colonne)

Organizzazione tabella doppia entrata

xtabs(frequenza~var1+var2) # frequenza sarà all'interno della tab, var1 righe, var2 colonne

Distribuzione di frequenza

xtabs(~x)

Fornire le percentuali di successo condizionate all’esplicativa var1

xtabs(successo~var1)/(xtabs(insuccesso~var1)+xtabs(successo~var1))

Odds ratio e rischio relativo

Calcolare l'odds ratio e il rischio relativo:

OR=tab[1,1]*tab[2,2]/(tab[1,2]*tab[2,1])
RR=tab[2,2]/tab[1,2] # gli anziani hanno RR volte la prob di malattia dei giovani

Si fornisca RR e OR del successo per var1 per uno studente con le altre esplicative pari alla media

new=data.frame(psi=c(0,1), tuce=c(mean(tuce), mean(tuce)), gpa=c(mean(gpa), mean(gpa)))
p=predict(mod, newdata=new, type="response")
RR=p[2]/p[1] # gli studenti sottoposti al nuovo metodo hanno una prob di migliorare i voti pari a RR volte quella degli studenti non sottoposti al nuovo metodo
OR=(p[2]/(1-p[2]))/(p[1]/(1-p[1]))

Fornire una previsione per la risposta quando l’esplicativa è pari a 31

ilogit(31) # binomial regression
log(31) # poisson regression

Fornire una previsione per la risposta per due nuovi ingressi aventi var1=max(var1) e var2="x1" e var2="x2"

new=data.frame(var1=max(var1), var2=c("x1","x2"))
predict(mod, new, type="response")
predict(mod, type="link") # stima il predittore lineare eta
predict(mod, type="response") # stima le probabilità theta (frequenze relative previste)

Stime del predittore lineare in corrispondenza dei valori di x osservati

eta=coef(mod)[1]+coef(mod)[2]*x
# per stimare la prob di successo per una singola unità statistica sostituisco ad x il valore dell'esplicativa di quella u.s
# stime della risposta media (prob successo)
theta=exp(eta)/(1+exp(eta)) # binomial
theta= exp(eta) # poisson

Confrontiamo i valori di Y osservati e le prob stimate

cbind(y, round(theta, 2)) # ovviamente la somma dei successi coincide con la somma delle prob previste

Riportare le probabilità stimate dal modello per ciascun pattern delle esplicative

cbind(dati, prev=round(fitted(mod), 2))

Riportare le frequenze assolute osservate della risposta e quelle previste dal modello

n=successi+insuccessi
freq.oss=successi
freq.attese=n*fitted(mod)
frequenze=round(cbind(freq.oss, freq.attese), 2)

Previsione intervallare

x=80
New=data.frame(Age=80) # devo specificare tutte le covariate del dataset (in questo ho solo age)
logit.80=predict(mod, newdata=new, type="link")
sigma.hat=summ$cov.unscaled # matrice di varianze e covarianze asintotica stimata
var.logit.hat=sigma.hat[1,1]+x^2*sigma.hat[2,2]+2*x*sigma.hat[1,2]
se.logit=sqrt(var.logit.hat)
IC.logit.80=logit.80+c(-1,1)*qnorm(0.975)*se.logit
IC.prob.80=exp(IC.logit.80)/(1+exp(IC.logit.80))

IC per beta

confint(mod) # valore d'indifferenza: 0

IC per OR

exp(confint(mod)) # valore d'indifferenza: 1

Costruire, per ciascuna esplicativa, un IC per l’odds ratio al 95% e trarne conclusioni

round(exp(cbind(coef(mod), confint(mod))), 3)

Associazione negativa e positiva fra un’esplicativa e la risposta

Associazione negativa: le stime puntuali dell’odds ratio sono minori di 1.
Associazione positiva: le stime puntuali dell’odds ratio sono maggiori di 1.

Trasformare un dataset con risposta binaria in binomiale e viceversa

Binario in binomiale:

patterns=prodotto del numero di modalità di ogni esplicativa (n° combinazioni possibili)
sum(pclass==1 & sex=='male' & age=='Giovane' & survived==0)
sum(pclass==1 & sex=='male' & age=='Giovane' & survived==1)
sum(pclass==1 & sex=='male' & age=='Anziano' & survived==0)
sum(pclass==1 & sex=='male' & age=='Anziano' & survived==1)
sum(pclass==1 & sex=='female' & age=='Giovane' & survived==0)
sum(pclass==1 & sex=='female' & age=='Giovane' & survived==1)
sum(pclass==1 & sex=='female' & age=='Anziano' & survived==0)
sum(pclass==1 & sex=='female' & age=='Anziano' & survived==1)