Statistica III - Teoria + Pratica

nulliLM0=as.numeric(logLik(mod.nullo))

LM=as.numeric(logLik(mod))

w.oss=-2*(LM0-LM)

p.value=1-pchisq(w.oss,diff.gdl)

#Il test di wald e LRT per confrontare due modelli che differiscono per un’esplicativa sono asintoticamente equivalenti,ma non sono lo stesso test

#e. P-value del test basato sulla differenza delle devianze:1-pchisq((mod$deviance-mod1$deviance),(mod$df.residual-mod1$df.residual))

#analogo ad anova(mod,mod1, test=’Chisq’) (non farti confondere dagli ***)

H0: i coefficienti aggiunti sono nulli

#TEST DI HOSMER-LEMESHOW (da usare in caso di dati sparsi)

X.HL= Σ (succ.oss-succ.att)^2/succ.att + Σ (insucc.oss-insucc.att)^2/insucc.att

pvalue=1-pchisq(X.HL,g-2) #H0: il modello è adeguato

#TEST PER LA VERIFICA DELLA SOVRADISPERSIONE

H0: assenza di sovradispersione

Toss=sum((data$risposta-(mean(data$risposta)))^2)/mean(data$risposta)

1-pchisq(Toss,(nrow(data)-1))

#SE LA DEVIANZA E' "TROPPO ELEVATA"

Prima di concludere che c’è

sovradispersione cerchiamo altre potenziali spiegazioni

#1. Dati sparsi: ripetizioni < 5 e/o successi stimati per l’i-esimo pattern < 1 (a occhio)data$risposta

#2. Correttezza modello (link, linearità, punti inusuali)par(mfrow=c(2,2))plot(mod)

#TRATTARE LA SOVRADISPERSIONE

#Si può:

1. stimare il parametro di dispersione e chiedere un summary del modello stimatophi=sum(residuals(mod,type="pearson")^2)/df.residual(mod)

#stima parametro di dispersione: X^2/gdl (nel modello poisson dovrebbe essere 1)summary(mod, dispersion=phi)

2. stimare il glm con family quasi-verosimiglianzasummary(glm(y~.,family=quasipoisson,data))

#le stime dei parametri non cambiano mentre vengono aggiustati gli SE e quindi la significatività osservata. Infatti gliSE vengono moltiplicati per sqrt(phi).NB: usando la quasi-verosimiglianza non posso più usare la devianza e tutte le qnt. che la coinvolgono (perde disignificato)

#FORNIRE IL PATTERN DELLE ESPLICATIVE CHE DÀ

Il maggior contributo alla devianza è dato da:

residuals(mod)^2

Le devianze dei pattern, se sommate, danno la devianza del modello.

Per valutare quale modello preferire tra mod e mod1 (modelli annidati), si può utilizzare il test:

anova(mod, mod1, test='Chisq')

Il test è basato sulla differenza di devianze e l'ipotesi nulla (H0) è che i coefficienti aggiunti nel modello più complesso siano nulli.

Per valutare i modelli basandosi su criteri, si può utilizzare il criterio AIC:

AIC(mod)

AIC(mod1)

Per calcolare il R^2 di Naglekerke tra modelli annidati, si può utilizzare la seguente formula:

n = sum(successi + insuccessi)
R2.NAG = (1 - exp((mod$dev - mod$null) / n)) / (1 - exp((-mod$null) / n))

Il valore di R2.NAG rappresenta la percentuale di spiegazione della varianza del modello corrente rispetto al modello saturo.

Per effettuare una selezione stepwise basata sull'AIC, si può utilizzare:

step(mod)

Le procedure backward e forward porteranno alla stessa conclusione.

Per fornire un test e un grafico per gli outlier, si può utilizzare la libreria "car" e i seguenti comandi:

outlierTest(mod)

plot(rstandard(mod))
abline(h = c(-2, 2))

Per fornire un test e un grafico per i punti di leva, si può utilizzare il seguente codice:

k = 4
soglia = 2 * (k + 1) / n
hatvalues(mod)[hatvalues(mod) > soglia]
plot(hatvalues(mod))

Disegna una linea orizzontale sul grafico al valore della soglia.

Trova l'indice del punto più influente utilizzando la distanza di Cook.

Mostra un grafico riassuntivo per outlier, punti di leva e punti influenti.

Calcola le misure di influenza per il modello.

Fornisce i residui del modello utilizzando il metodo di devianza.

Fornisce i residui del modello utilizzando il metodo di Pearson.

Fornisce i residui del modello in termini di risposta: frequenza relativa osservata - probabilità stimata.

Fornisce un tipo di residui sotto prodotto della procedura IRLS, non servono per la diagnostica.

Calcola la matrice H di leva del modello, che è una funzione sia delle variabili esplicative che della risposta.

La somma dei valori della matrice H è uguale a k+1, dove k è il numero di variabili esplicative.

Calcola i residui standardizzati e studentizzati del modello utilizzando le stesse funzioni degli lm.

residui standardizzati (di devianza o di pearson)

#SI FORNISCA UN GRAFICO UTILE A VALUTARE L’APPROPRIATEZZA DELLA FUNZIONE DI VARIANZA

aspettativa: dispersione costante

plot(fitted(mod),rstandard(mod))

#DIAGNOSTICA GRAFICA DI BASE FORNITA DA R PER IL MODELLO STIMATO

Non funziona nel caso di dati sparsi

par(mfrow=c(2,2))

plot(mod)

#GRAFICO 1: RESIDUALS VS FITTED. Valuto l’evidenza di non linearità (eta.hat vs residui di devianza (eD))

#GRAFICO 2: NORMAL Q-Q. Valuto la normalità asintotica dei residui. Qualche osservazione da indagare meglio?

#GRAFICO 3: SCALE-LOCATION. Valuto problemi nella funzione di varianza (eta.hat vs sqrt(eDS))

#GRAFICO 4: RESIDUALS VS LEVERAGE. Valuto l’anomalia di singole osservazioni (h.ii vs ePS)

a. i residui usati sono quelli di devianza

b. i fitted values sono in termini di predittore lineare g(mui) e non di valori previsti per la risposta (mui)

c. un leggero allontanamento dalla bisettrice nel QQ-plot e' frequente

Se valuto problemi diagnostici,

segmenti ottenuti non sono paralleli potrebbe essere necessario inserire nel mod l'interazione tra le 2 esplicative

#PLOT DELLE FREQUENZE RELATIVE AL VARIARE DI VAR1

plot(y.i~var1)

#PLOT DEL MODELLO STIMATO SOVRAIMPOSTO ALLE PROPORZIONI OSSERVATE

plot(freq.relativa.risposta~var1)

lines(var1,fitted(mod))

#PLOT DATI CON SOVRAIMPOSTA LA CURVA DEL MODELLO STIMATO

plot(x,y)

a=seq(0,90,0.5)

b=coef(mod)[1]+coef(mod)[2]*a #a moltiplica il coeff di x

lines(a,ilogit(b), type='l') #Se binomial

lines(a,exp(b), type='l') #se poisson

#Si disegni un grafico riportante la curva stimata al variare dei valori dell'esplicativa var1 (caso A) per uno studente che non è stato esposto al nuovo metodo di insegnamento e ha valore della variabile tuce pari alla sua media e (caso B) per uno studente che è stato esposto al nuovo metodo di insegnamento e ha valore della variabile tuce pari alla sua media. Il range per l'esplicativa deve essere (2,4).

plot(risposta~var1, xlim=c(2,4),ylim=c(0,1),pch=20)

x=seq(2,4,by=.1)

casoA=coef(mod)[1]+coef(mod)[3]*mean(tuce)+coef(mod)[4]*x

casoB=coef(mod)[1]+coef(mod)[2]+coef(mod)[3]*mean(tuce)+coef(mod)[4]*x

lines(x,

ilogit(casoA),type="l",lty=2)
lines(x, ilogit(casoB),type="l",lty=3)
legend(2,0.9,lty=c(2,3),legend=c("psi=0","psi=1"))
#Nel modello saturo la devianza, i residui di devianza e i gdl residui risultano tutti nulli
#la devianza residua di un lm coincide con la devianza di un glm
#Se un modello non mi convince (devianza ancora alta rispetto ai gdl) posso sempre provare tramite trasformate delle esplicative
#Numero di iterazioni del Fisher-Scoring: 4 dopo l'inizializzazione l'algortimo numerico di stima necessita di sole 3 iterazioni
#ZERO-INFLATED POISSON REGRESSION MODEL: trattamento della sovradispersione