vuoi
o PayPal
tutte le volte che vuoi
DECIDERE TRA DUE MODELLI CON 2 TEST DIVERSI
nova(mod, mod2, test='Chisq') #confronto con il modello precedenteà#il p-value è quasi zero rifiutiamo il modelloprecedente in favore di quello nuovof.obs <- (dev-dev2)/(dev2/gdl2) #alternativa all'ANOVA ma è non precisa1-pf(f.obs,1,gdl2)AIC(mod2) #per il criterio AIC scegliamo il modello che ha AIC minoreAIC(mod3)
DECIDERE SE UNA VARIABILE RISPOSTA È NORMALMENTE DISTRIBUITA
shapiro.test (Y) #testqqnorm (Y) #graficamentecurve(dnorm(x, mean(Y), sd(Y), add=T)plot (density(Y)) #motivazione teorica
DECIDERE SE UNA VARIABILE RISPOSTA È DISTRIBUITA COME UNA POISSON
#distribuzione empiricatab <- xtabs(~Y, data=data)ascisse <- as.numeric(names(tab))plot(ascisse, tab)points(ascisse, tab) #dist. empirica = dist. teorica allora penso che Poisson non vada bene#distribuzione teoricateorica <- rpois(250, mean(Y))tabT<- xtabs( ~ teorica)ascisseT <- as.numeric(names(tabT))plot(ascisseT, tabT)points(ascisseT,
```htmlT)ks.test(jitter(tab), "ppois", mean(Y)) #test
library(car)
qqPlot(tab, distribution="pois", lambda=mean(Y)) #motivazione teorica
CROSS PRODUCT RATIO E RISCHIO RELATIVO
tab <- xtabs(Y~X +X , data=data)
1 2
OR=tab[1,1]*tab[2,2]/(tab[1,2]*tab[2,1])
RR=tab[2,2]/tab[1,2]
SI FORNISCA MISURA RR DEL SUCCESSO (Y=1) PER UNA VARIABILE 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]
OR <- (p[2]/(1-p[2]))/(p[1]/(1-p[1])) à
COMMENTO I PLOT DI DIAGNOSTICA RESTITUITI DA R restituisce 4 grafici
par(mfrow=c(2,2)) #per vedere tutti e 4 i grafici insieme
plot(mod)
Grafico 1: per verificare la linearità del modelloà eta.hat VS residui di devianza
Grafico 2: QQplot per verificare la presenza di eventuali outlier
Grafico 3: serve a verificare la validità della funzione di varianza che
```stiamo utilizzando che risulta appunto valida se abbiamo un andamento costante
Grafico 4: punti estremi verso destra indicano punti di leva, mentre punti estremi verso l'alto o il basso indicano potenziali outlier h.
ii VS residui di pearson standardizzati
Se valuto problemi diagnostici, come prima cosa cerco di risolvere la linearità aggiungendo le esplicative al quadrato/cubo/log in base al trend del primo grafico, es: I(var^2)
Così facendo potrei risolvere altri problemi diagnostici.
Se ad esempio riesco a correggere la linearità trovando un modello che fitta molto bene, posso non preoccuparmi se la funzione di varianza non è buona. Elimino osservazioni come ultima spiaggia.
DIAGNOSTICA PER FORMA STRUTTURALE: Linearità. Plot eta.hat (predittore lineare stimato) e e^DS (residui di devianza standardizzati)
Mi aspetto linearità: plot(mod$lin, rstandard(mod)) #se non lineare aggiungo nel modello un coeff: predict(mod)^2 e valuto se
è significativoà Link. Plot di eta.hat vs zeta.hat (stima pseudo-risposte)
Mi aspetto un andamento rettilineo, altrimenti probabilmente il link è inopportuno
library(faraway)
eta.hat <- predict(mod)
mu.hat <- ilogit(eta.hat) #Binomial regression
mu.hat <- exp(eta.hat) #Poisson regression
derivata <- 1/(mu.hat*(1-mu.hat)) #Binomial regression
derivata <- 1/(mu.hat) #Poisson regression
zeta.hat <- eta.hat + (Y-mu.hat)*derivata
plot(zeta.hat ~ eta.hat)à Funzione di varianza. mu.hat (risposta media stimata) vs residui standardizzati (didevianza o di pearson)
Mi aspetto dispersione costante
par(mfrow=c(1,1))
mu.hat <- fitted(mod)
plot(rstandars(mod) ~ mu.hat)
plot(residuals(mod, type=”responde”) ~ mu.hat)DATO LO Z-VALUE CALCOLO IL P-VALUE
z.value <- summary(mod)$coef[1,3]
2*(1-pnorm(abs(z.value)))DATO IL P-VALUE CALCOLO LO Z-VALUE
p.value <- summary(mod)$coef[1,4]
qnorm(1-p.value/2)*(-1)GRAFICO CURVA STIMATA(PLOT) SOVRAIMPOSTA AI PUNTI
OSSERVATI..range(30,85)par(mfrow=c(1,1))plot(Y~ X +X , xlim=c(30,85), ylim=c(0,1))1 2grid <- seq(30,85,by=1)eta.logit <- coef(mod)[1]+coef(mod)[2]*gridmu.logit <- ilogit(eta.logit) #Binomialmu.logit <- exp(eta.logit) #Poissonlines(grid, mu.logit, type="l")STIMA DEL PARAMETRO DI DISPERSIONEphi.hat <- sum(residuals(mod, type=”pearson”)^2)/df.residuals(mod)summary(mod, dispersion=phi.hat)TEST PER LA VERIFICA DELLA SOVRADISPERSIONEH : assenza di sovradispersione0n <- nrow(dataset)test <- sum(Y-(mean(Y))^2)/mean(Y)1-pchisq(test , (n-1))TRATTARE LA SOVRADISPERSIONE1. stimare il parametro di dispersione e chiedere un summary del modello stimatophi <- sum(residuals(mod, type="pearson")^2)/df.residual(mod)summary(mod, dispersion=phi)2. stimare il glm con family=quasi-verosimiglianzamod <- glm(Y~., family=quasipoisson, data=data)summary(mod)#le stime dei parametri non cambiano mentre vengono aggiustati gli SE e quindi la significativitàosservata.
Infatti gli SE vengono moltiplicati per sqrt(phi).NB: usando la quasi-verosimiglianza non posso più usare la devianza e tutte le quantità che la coinvolgono poiché perde di significato. PREVISIONE PER LA Y QUANDO X =a , X =b , X =mean(X )<new - data.frame(X =a , X =b , X =mean(X ))
<predict(mod, new, type="response")
IC per beta:
<confint(mod) #valore d’indifferenza: 0
IC per OR:
<exp(confint(mod)) #valore d’indifferenza: 1
SELEZIONE BACKWARD SIA CON CRITERION-BASED CHE CON TEST-BASED
<mod.back - step(mod, method="backward") #criterion-based
<summary(mod.back)
<summary(mod)
<mod2 - update(mod,.~. -X) #test-nased
<summary(mod2 )
VALUTARE OUTLIER
<library(car)
<res - rstandard(mod)
<res.[abs(res)>2] #approssimativo
<outlierTest(mod) #aggiustato secondo Bonferroni
<plot(rstandard(mod))
<abline(h=c(-2,2))
VALUTARE LEVERAGE
<libary(car)
<soglia - 2*(k+1)/n
<h.ii -<code>hatvalues(mod)[hatvalues(mod)>soglia]</code>
<code>plot(hatvalues(mod))</code>
<code>abline(h=soglia)</code>
<code>influence.measures(mod)</code>
<code>VALUTARE OSSERVAZIONE PIU' INFLUENTE</code>
<code>library(car)</code>
<code>which.max(cooks.distance(mod))</code>
<code>influencePlot(mod)</code>
<code>influenceIndexPlot(mod)</code>
<code>TRASFORMARE UN DATASET CON RISPOSTA BINARIA IN BINOMIALE E VICEVERSA</code>
<code>Binario in binomiale: (1 riga di codice)</code>
<code>patterns=prodotto del numero di modalità di ogni esplicativa (n° combinazioni possibili)</code>
<code>sum(pclass==1 & sex=='male' & age=='Giovane' & survived==0)</code>
<code>sum(pclass==1 & sex=='male' & age=='Giovane' & survived==1)</code>
<code>sum(pclass==1 & sex=='male' & age=='Anziano' & survived==0)</code>
<code>sum(pclass==1 & sex=='male' & age=='Anziano' & survived==1)</code>
<code>sum(pclass==1 & sex=='female' & age=='Giovane' & survived==0)</code>
<code>sum(pclass==1 & sex=='female' & age=='Giovane' & survived==1)</code>
<code>sum(pclass==1
<p>& sex=='female' & age=='Anziano' & survived==0)sum(pclass==1 & sex=='female' & age=='Anziano' & survived==1)sum(pclass==2 & sex=='male' & age=='Giovane' & survived==0)sum(pclass==2 & sex=='male' & age=='Giovane' & survived==1)sum(pclass==2 & sex=='male' & age=='Anziano' & survived==0)sum(pclass==2 & sex=='male' & age=='Anziano' & survived==1)sum(pclass==2 & sex=='female' & age=='Giovane' & survived==0)sum(pclass==2 & sex=='female' & age=='Giovane' & survived==1)sum(pclass==2 & sex=='female' & age=='Anziano' & survived==0)sum(pclass==2 & sex=='female' & age=='Anziano' & survived==1)sum(pclass==3 & sex=='male' & age=='Giovane' &
Per specificare la variabile risposta:
- Si fornisce la proporzione di successi (yi=x.i/n.i) come risposta e, tra gli argomenti della funzione glm, si specificano il numero di repliche per ogni valore delle esplicative in weights (n.i ≠ n)
- Si fornisce la risposta come array costituito da due colonne: la prima contiene il numero di successi (yitilde) e la seconda il numero degli insuccessi (ni-yitilde)
mod <- glm(y.i~X +X +…+X , family=binomial, weights =n.i, data)mod <- glm(cbind(successo,ni-successo)~., family=binomial, data)Il binary model si implementa come il binomial.
È possibile effettuare un test F basato sulla differenza delle devianze e un test Hosmer-Lemeshow.
Non è possibile effettuare un test per la bontà del modello, residui di devianza e diagnostica poiché i dati sono sparsi.
Quindi, è necessario trasformare il dataset in modo da avere per ogni riga un pattern (binomial regression).
L'equazione del modello stimato è: mu=exp(.) / 1+exp(.).
Per stimare un GLM, si fa ricorso al link canonico.
<p>(POISSON)mod1 <- glm(y~., family=poisson, data)
mod2 <- glm(y~offset(log(var1)) + var2 + … + varj, family=poisson, data)
Inserisco una variabile offset in taglia logaritmica se devo ovviare alle differenze di sizeà dell’esplicativa: n° esposizioni molto spesso me ne accorgo perché tutti i B sono mooolto piccoli.
Se non posso perché non ho la variabile n° di esposizioni ma noto valori decisamente diversi nel dataset in qualche esplicativa (non risposta!!) (es: 0.5 per una u.s. e 4000 per un’altra), provo a trasformarle es: log mod=glm(y~var.qual*var1.quant+var.qual*v)
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo
