# Codici in R per la regressione logistica # Importo i dati dei pozzi ##invlogit<- function (x) ##{ ## 1/(1 + exp(-x)) ##} setwd("/Users/mariagraziapittau/Library/CloudStorage/Dropbox-Universita'DegliStudiDiRomaLaSapienza/grazia.pittau@uniroma1.it/Lezioni/Statistica_Aziendale/R//Logistico") library(arm) wells <- read.table ("wells.dat", header=TRUE) wells[1:5,] attach (wells) # Regressione logistica relativa alla probabilita' # di cambiare pozzo per il rifornimento dell'acqua # in funzione della SOLA distanza dal pozzo sicuro piu' vicino fit.1 <- glm (switch ~ dist, family=binomial(link="logit")) display (fit.1,5) # Istogramma della distanza hist (dist, breaks=seq(0,10+max(dist[!is.na(dist)]),10), freq=T, xlab="Distanza (in metri) dal pozzo sicuro più vicino", ylab="", main="", mgp=c(2,.5,0)) # Ridefinisco la distanza esprimendola in 100 metri (distanza a base 100) dist100 <- dist/100 fit.2 <- glm (switch ~ dist100, family=binomial(link="logit")) display (fit.2) # Rappresentazione del modello stimato # Funzione di jitter per dati binari # dati che sono del tipo 0 e 1 jitter.binary <- function(a, jitt=.05){ a + (1-2*a)*runif(length(a),0,jitt) } ##La funzione beta.hat per estrarre i coefficienti beta.hat <- function(fit){ summ <- summary (fit, correlation=T) beta.hat <- summ$coef[,1] } ## Da notare: # 1. la funzione jitter definita per dati che assumono solo valori del # tipo 0 e 1 # La funzione invlogit plot(c(0,max(dist100, na.rm=T)*1.02), c(0,1), xlab="Distanza (in 100 metri) dal pozzo sicuro più vicino", ylab="Pr (cambiare)", type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0)) points (dist100, jitter.binary(switch), pch=20, cex=.1) curve(invlogit(coef(fit.2)[1]+coef(fit.2)[2]*x), lwd=1, add=T) ##La prob di cambiare pozzo scende da 0.6 per chi abita vicino a # 0.2 per le famiglie che vivono lontane da un pozzo sicuro piu' di 300 metri ##Valore del predittore lineare nel punto medio pred.lin<-beta.hat(fit.2)[1]+beta.hat(fit.2)[2]*mean(dist100) ##calcolo la pendenza della curva in un intorno infinitesimale del punto medio d.dx<-beta.hat(fit.2)[2] * exp(pred.lin) /(1+exp(pred.lin))^2 reg.4<-beta.hat(fit.2)[2]/4 ## buona approssimazione perche' la curva nel punto medio pari a 0.4 passa vicino al 50% e intorno a questo valore si trovano la maggiorparte dei dati (siamo intorno al punto di massima pendenza)---se infatti considero 250 metri invece del valor medio l'approssimazione meggiora notevolmente. ##Aggiunta di un predittore: # Livello di concentrazione di arsenico nel pozzo fit.3 <- glm(switch ~ dist100 + arsenic, family=binomial(link="logit")) display (fit.3) (beta.hat(fit.3)[2])/4*100 ##22% (beta.hat(fit.3)[3])/4*100 ##12% ## la distanza ha un'importanza superiore alla quantita' di arsenico--risultato inverosimile ##CALCOLO LE SD sd(dist100) sd(arsenic) ## La distanza e' molto meno variabile del livello di arsenico ## Il problema e' che DEVO standardizzare ##calcolo i coefficienti corrispondenti a una volta la deviazione standard ##La variabilita' della distanza e' molto piu' bassa di quella dell'arsenico beta.hat(fit.3)[2]*sd(dist100) beta.hat(fit.3)[3]*sd(arsenic) (beta.hat(fit.3)[2]*sd(dist100))/4*100 ##8% (beta.hat(fit.3)[3]*sd(arsenic))/4*100 ##13% ##L'effetto del veleno e' molto piu' importante dell'effetto della distanza ##SOLUZIONE: STANDARDIZZAZIONE fit.a <- glm(switch ~ scale(dist100) + scale(arsenic), family=binomial(link="logit")) display (fit.a) ##Rappresentazione grafica (con il modello senza riscalare per convenienza visiva) # In funzione della distanza plot(c(0,max(dist,na.rm=T)*1.02), c(0,1), xlab="Distanza (in metri) dal pozzo sicuro più vicino", ylab="Pr (cambiare)", type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0)) points (dist, jitter.binary(switch), pch=20, cex=.1) curve (invlogit(beta.hat(fit.3)[1]+beta.hat(fit.3)[2]*x/100+beta.hat(fit.3)[3]*.52), lwd=2.5, add=T, lty=1) curve (invlogit(beta.hat(fit.3)[1]+beta.hat(fit.3)[2]*x/100+beta.hat(fit.3)[3]*2.50), lwd=2.5, add=T, lty=2, col="red") curve (invlogit(beta.hat(fit.3)[1]+beta.hat(fit.3)[2]*x/100+beta.hat(fit.3)[3]*6.50), lwd=5.5, add=T, lty=3, col="pink") text (50, .65, "se Arsenico = 2.50", adj=0, cex=1.8) text (50, .45, "se Arsenico = 0.52", adj=0, cex=1.8, col="red") text (50, .55, "se Arsenico = 6.5", adj=0, cex=1.8, col="pink") #In funzione del livello di arsenico plot(c(0,max(arsenic,na.rm=T)*1.02), c(0,1), xlab="Concentrazione di arsenico nell'acqua dei pozzi", ylab="Pr (cambiare)", type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0)) points (arsenic, jitter.binary(switch), pch=20, cex=.1) curve (invlogit(beta.hat(fit.3)[1]+beta.hat(fit.3)[2]*0+beta.hat(fit.3)[3]*x), from=0.5, lwd=2.5, add=T, col="red") curve (invlogit(beta.hat(fit.3)[1]+beta.hat(fit.3)[2]*0.5+beta.hat(fit.3)[3]*x), from=0.5, lwd=2.5, add=T, col="black") curve (invlogit(beta.hat(fit.3)[1]+beta.hat(fit.3)[2]*1+beta.hat(fit.3)[3]*x), from=0.5, lwd=2.5, add=T, col="blue") text (6, .7, "se distanza = 0 metri", adj=0, cex=1.8, col="red") text (6, .65, "se distanza = 50 metri", adj=0, cex=1.8) text (6, .6, "se distanza = 100 metri", adj=0, cex=1.8, col="blue") #Modello con interazione fit.4 <- glm (switch ~ scale(dist100) + scale(arsenic) + scale(dist100):scale(arsenic), family=binomial(link="logit")) display(fit.4) ##Valuto tutto in termini di punti medi..il valore dell'arsenico pari a 0 non ha senso perche' il valore soglia e' 0.5 ##Il valore del predittore lineare nei punti medi pred.lin <-beta.hat(fit.4)[1] - beta.hat(fit.4)[2] * mean(dist100) + beta.hat(fit.4)[3] * mean(arsenic) + beta.hat(fit.4)[4]* mean(dist100)*mean(arsenic) invlogit(pred.lin) ##-.15 -.58*.48 + .56*1.66-.18*.48*1.66 pred.lina <-beta.hat(fit.a)[1] - beta.hat(fit.a)[2] * mean(dist100) + beta.hat(fit.a)[3] * mean(arsenic) invlogit(pred.lina) ##Rappresentazione grafica del modello con interazione (valori grezzi per semplicita') ##In funzione dell'arsenico plot(c(0,max(dist,na.rm=T)*1.02), c(0,1), xlab="Distanza (in metri) dal pozzo sicuro più vicino", ylab="Pr (cambiare)", type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0)) points (dist, jitter.binary(switch), pch=20, cex=.1) curve (invlogit (cbind(1,x/100, 9, 9*x/100) %*% coef(fit.4)), add=TRUE, lwd=3) curve (invlogit (cbind(1,x/100,.52,0.52*x/100) %*% coef(fit.4)), add=TRUE, col="red", lwd=3) text (75, .750, "se Arsenico = 9.0", adj=0, cex=1.8) text (50, .49, "se Arsenico = 0.52", adj=0, cex=1.8, col="red") #In funzione della distanza plot(c(0,max(arsenic,na.rm=T)*1.02), c(0,1), xlab="Livello di concentrazione di arsenico", ylab="Pr (cambiare)", type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0)) points (arsenic, jitter.binary(switch), pch=20, cex=.1) curve (invlogit(beta.hat(fit.4)[1]+beta.hat(fit.4)[2]*0+beta.hat(fit.4)[3]*x+beta.hat(fit.4)[4]*1*x), from=0.5, lwd=3.5, add=T, col="red") curve (invlogit(beta.hat(fit.4)[1]+beta.hat(fit.4)[2]*0.5+beta.hat(fit.4)[3]*x+beta.hat(fit.4)[4]*0.5*x), from=0.5, lwd=3.5, add=T, col="blue") curve (invlogit(beta.hat(fit.4)[1]+beta.hat(fit.4)[2]*0.9+beta.hat(fit.4)[3]*x+beta.hat(fit.4)[4]*0.9*x), from=0.5, lwd=.5, add=T, col="pink") text (6, .8, "se distanza = 1", adj=0, cex=.1.8, col="red") text (6, .6, "se distanza = 50", adj=0, cex=1.8, col="blue") text (6, .9, "se distanza = 90", adj=0, cex=1.8, col="pink") #Aggiunta deipredittori sociali ##educ===anni di istruzione dell’utilizzatore del pozzo. educ4 <- educ/4 ##rappresenta la differenza attesa conseguente all’aggiunta di quattro anni di istruzione ##assoc partecipazione o meno della famiglia a qualche organizzazione all'interno della comunita' fit.5 <- glm (switch ~ rescale(dist100) + rescale(arsenic) + rescale(dist100):rescale(arsenic) + assoc + rescale(educ4) + assoc:rescale(educ4), family=binomial(link="logit")) display (fit.5) # L'associazione non ha senso (il segno e' irragionevole) # e non risulta neanche significativa # la levo dal modello fit.6 <- glm (switch ~ scale(dist100) + scale(arsenic) + scale(dist100):scale(arsenic) + scale(educ4), family=binomial(link="logit")) display (fit.6) # OSS. La variabile educ ha un effetto molto alto, allora # come di consuetudine proviamo a considerare le interazioni di questa variabile # e altri predittori #Introduco le interazioni ##educ==anni di istruzione dell'utilizzatore del pozzo #Il nuovo modello con l'interazione tra educ e altri predittori fit.7 <- glm (switch ~ scale(dist100) + scale(arsenic) + scale(educ4) + scale(dist100):scale(arsenic) + scale(dist100):scale(educ4) + scale(arsenic):scale(educ4), family=binomial(link="logit")) display (fit.7) pred.7<-fit.7$fitted.values # 0.35 ma il coeff della distanza era negativo # un livello di educazione superiore riduce l'effetto negativo della distanza--indipendentemente dall'arsenico # 0.08 un aumento nel livello di educazione aumenta l'effetto positivo dell'arsenico ##ANALISI DEI RESIDUI #Calcolo dei residui ##Valori stimati pred.6 <- fit.6$fitted.values ##valori osservati y <- switch ##residui res <- y-pred.6 #Grafico dei residui grezzi plot(c(0,1), c(-1,1), xlab="Prob (cambiare ) stimata", ylab="Osservati - stimati", type="n", main="Grafico dei residui", mgp=c(2,.5,0)) abline (0,0, col="gray", lwd=.5) points (pred.6, y-pred.6, pch=20, cex=.2) par(mfrow=c(2,1)) #Grafico dei residui in funzione delle probabilita' stimate binnedplot (pred.6, y-pred.6, nclass=40, xlab="Prob (cambiare ) stimata", ylab="Residui medi", main="Grafico dei residui ''binned''", mgp=c(2,.5,0)) #Grafico dei residui in funzione di una variabile di input: l'arsenico binnedplot (arsenic, y-pred.6, nclass=80, xlab="Livello di arsenico", ylab="Residui medi", main="Grafico dei residui ''binned''", mgp=c(2,.5,0)) # Si nota un comportamento anomalo: residui positivi nella parte centrale e residui negativi per valori # elevati dell'arsenico ## 2 valori molto bassi per livello di arsenico circa 0.52 ##il valore previsto dello switching del modello e' 49% ma solo il 32% cambiano pozzo #Possibile soluzione: trasformazione dell'arsenico su scala logaritmica ##Modello nullo==modello con sola intercetta fit.0<- glm (switch ~ 1, family=binomial(link="logit")) display (fit.0) ##I valori stimti dal modello nullo pred.0<-fit.0$fitted.values table(pred.0) round(unique(pred.0), 3) ## la prob e' maggiore di 0.5 quindi il modello nullo assegna 1 a tutti i valori ##percentuale di 1 nei dati p.hat<-table(y)[2]/length(y) ##Tasso di errore nel modello nullo wr0<-1-p.hat table(y) ##proporzione di zero nei dati table(y)[1]/length(y) ##proporzione di 1s nei dati table(y)[2]/length(y) ##il 58% degli intervistati dichiara di essere favorevole a cambiare pozzo ##tasso di errore tax.error.null <- (length(y)-sum(y[y==1]))/length(y) error.rate.null <- mean(round(abs(y-mean(pred.6)))) ##Stimo il modello (fit.9) ##TASSO DI ERRORE==(no1+n10)/n mean(pred.6 > 0.5 & y==0) + mean(pred.6 <0.5 & y==1) ##ovvero scritto in un altro modo tax.error <- mean((pred.6 > 0.5 & y==0) | (pred.6 < 0.5 & y==1)) error.rate <- mean(round(abs(y-pred.6))) #media (y-y.estim) ## Calcolo della tabella a doppia entrata con valori osservati e valori stimati table(round(pred.6)) table(y.tilde=round(pred.6), y.oss=y) #DEVIANZA fit.1 <- glm (switch ~ 1, family=binomial(link="logit")) display (fit.1) #deviance = 4076.2 #Aggiungo un rumore newpred <- rnorm (length(switch),0,1) fit.1a <- glm (switch ~ newpred, family=binomial(link="logit")) display (fit.1a) #deviance = 4076.1 #Aggiungo un predittore che non sia un rumore casuale fit.1b <- glm (switch ~ dist100, family=binomial(link="logit")) display (fit.1b) # deviance = 3930.7 e' diminuita piu' di 1 ##Logit o probit fit.logit <- glm (switch ~ dist100 + arsenic, family=binomial(link="logit")) display (fit.logit) fit.probit <- glm (switch ~ dist100 + arsenic, family=binomial(link="probit")) display (fit.probit) ##Coeff probit trasformati in logit e interpretati con la regola del beta/4 ## probit*1.6=logit summary(fit.probit)$coef[2,1]*1.6 summary(fit.logit)$coef[2,1] ##Distribuzione normale e logistica curve(exp(x)/(1+exp(x))^2, -6, +6, ylab="", lty=1) curve(dnorm(x, 0, 1.6), add=TRUE, lty=2,lwd=2, -6, +6) legend(3, 0.15, c("Logistica", "Gaussiana"), lty=c(1,2))