########################################################################## # Laboratorio di "Metodi Statistici avanzati" con R # Prof.ssa Lea Petrella # Facolta' di Economia # Dipartimento Metodi e Modelli per l'Economia il Territorio e la Finanza # Sapienza Universita' di Roma # a.a. 2015-2016 # alberto.diiorio@uniroma1.it ########################################################################## ########################################################################## # Lezione 5: MODELLO di REGRESSIONE LINEARE MULTIPLA ########################################################################## # # Indice: # # * richiami teorici # * dati simulati dal "vero modello" # * analisi dei residui # * selezione delle variabili # * multicollinearit?† # * regressione multipla + variabili dummy # * dati reali (reddito) # * dati reali (circonferenza cranica) # ########################################################################## #---------------------------------------------------------- # * richiami teorici #---------------------------------------------------------- # p variabili esplicative (ovvero p-1 regressori effettivi + la costante unitaria) # epsilon e' il termine di errore...(aleatorio) ---> Y aleatoria # Y[i]=b_1+b_2X_2[i]...+b_kX_k[i]+epsilon[i] # il modello in forma matriciale: # Y=Xb+epsilon # Y vettore n x 1 # X matrice n x p # b vettore p x 1 # epsilon vettore n x 1 # assunzioni di base: # E(epsilon[i])=0 # V(epsilon[i])=sigma^2 # COV(epsilon[i],epsilon[j])=0 # Le ipotesi sul termine di errore possono essere sintetizzate # E(epsilon)=0 (vettore n x 1 di tutti zeri) # V(epsilon)=sigma^2*I (I matrice identit?† n x n) # ---> # E(Y|X)=Xb # V(Y|X)=sigma^2*I # Numero totale dei parametri ?® p+1 (parametri di regressione + varianza dell'errore) # generalizzate le assunzioni del modello lineare semplice # Matrice X di rango pieno: # nessuna delle colonne `e combinazione lineare delle altre # ---> X'X invertibile (det(X'X) diverso da zero) #---------------------------------------------------------- # * dati simulati dal "vero modello" #---------------------------------------------------------- n=100 p=3 X=matrix(0,nrow=n,ncol=p) b1=7 b2=2 b3=-3 b=c(b1,b2,b3) sigma2=100 set.seed(4065) X[,1]=1 X[,2]=round(runif(n,50,100),1) X[,3]=round(runif(n,-20,20),1) y=round(rnorm(n,X%*%b,sqrt(sigma2)),1) y dati=data.frame(X[,-1],y) dati round(cor(dati),2) cor.mat=cor(dati) library(corrplot) corrplot(cor.mat,method="circle") # Testiamo la significativita' del coefficiente di correlazione library(psych) corr.test(dati) # Verifiachiamo anche mediante l'ispezione grafica plot(dati$X1,dati$X2,xlab=expression(X[1]),ylab=expression(X[2]),main="Scatterplot tra i due regressori") plot(dati$X1,dati$y,xlab=expression(X[1]),ylab="y") plot(dati$X2,dati$y,xlab=expression(X[1]),ylab="y") # In un colpo solo faccio tutti gli scatter plot ma attenzione, # il comando ha senso se le variabili sono tutte quantitative! pairs(dati) plot(dati) #e' equivalente! # verifichiamo il rango della matrice sia pieno # nel nostro caso che sia uguale a 3 library(Matrix) rankMatrix(X)[[1]] # variabili esplicative indipendenti library(scatterplot3d) library(class) library(car) scatter3d(y~X[,2]+X[,3],surface=F,xlab="X2",ylab="Y",zlab="X3") # Richiede Pacchetto rgl # Stime dei minimi quadrati # hat.b=(X'X)^(-1) X'Y ---> hat.Y=X*hat.b b.hat=solve(t(X)%*%X)%*%t(X)%*%y b.hat # Gauss-Markov: hat.b stimatore BLUE # valori attesi ovvero stimati dall'iperpiano di regressione y.hat=X%*%b.hat #Restituisce una matrice di una colonna e 100 righe dim(y.hat) y.hat=c(y.hat) #Vettore di 100 componenti # residui residui=y-y.hat # un residuo elevato in valore assoluto indica che l'unita' # statistica e' descritta male dal modello # scomposizione della devianza totale SST=var(y)*(n-1) SSR=var(y.hat)*(n-1) SSE=var(residui)*(n-1) SST SSR+SSE # stima di sigma2 sigma2.hat=SSE/(n-p) sigma2.hat # standard error per sigma^2 sqrt(sigma2.hat) # coefficiente di determinazione e corrispondente versione aggiustata R2=1-SSE/SST R2 R2corretto=1-SSE/SST*(n-1)/(n-p) R2corretto # stima della matrice di varianze e covarianze delle stime b.hat vcov.mat.b.hat=sigma2.hat*solve(t(X)%*%X) vcov.mat.b.hat # stime delle varianze (e deviazioni standard) delle stime b.hat var.b.hat=diag(vcov.mat.b.hat) # stima della varianza se.b.hat=sqrt(var.b.hat) # stima della deviazione standard = STANDARD ERROR! se.b.hat # NOTA: finora nessuna ipotesi sulla distribuzione di epsilon! # aggiungiamo l'ipotesi di normalita' degli errori # ---> stesse stime dei minimi quadrati per a e b # test di significativita' t.oss.b=b.hat/se.b.hat t.oss.b # calcoliamo i p-value (occhio al valore assoluto...) p.value.b=2*pt(abs(t.oss.b),n-p,lower.tail=F) p.value.b # ANOVA: test F sui coefficienti (simultaneo) # H0: b2 = ... = bp = 0 # H1: altrimenti F.oss=(SSR*(n-p))/(SSE*(p-1)) F.oss pf(F.oss,df1=(p-1),df2=(n-p),lower.tail=F) # ...tutto piu' semplice con la funzione lm() di R mod=lm(y~X1+X2,data=dati) summary(mod) # Rappresentiamo il piano di regressione scatter3d(y~dati$X1+dati$X2,surface=T,xlab="X2",ylab="y",zlab="X3") # Intervalli di confidenza per i parametri di regressione confint(mod) # ora vediamo l'intervallo per la varianza dell'errore sigma2.int95=c((n-p)*sigma2.hat/qchisq(p=0.975,df=n-p),(n-p)*sigma2.hat/qchisq(p=0.025,df=n-p)) sigma2.int95 #---------------------------------------------------------- # * analisi dei residui #---------------------------------------------------------- predetti=mod$fitted.values residui=mod$residuals plot(residui) abline(h=0) # Verifichiamo che la media sia nulla con un opportuno test t.test(residui) # plottiamo i residui plot(predetti,residui,xlab=expression(hat(y)),ylab="residuals") abline(h=0) par(mfrow=c(2,2)) plot(mod) #test omoschedasticita'-eteroschedasticita' ncvTest(mod) library(lmtest) bptest(mod) # Per l'incorrelazione si puo' costruire il # correlogramma acf(residui) # il test di Ljung-Box Box.test(residui,lag=1,type="Ljung-Box") #Verifico per lag di ordine 1 Box.test(residui,lag=2,type="Ljung-Box") #Verifico per lag di ordine 2 # e test di Durbin-Watson dwtest(mod) # Verifichiamo la normalita' (Analisi Qualitativa) hist(residui,prob=T,xlim=c(-50,50)) curve(dnorm(x,0,sd(residui)),add=T,col="red",lwd=2) # Costruiamo il qqplot per i residui standardizzati (Analisi Qualitativa) qqnorm(scale(residui)) abline(0,1) # Test di normalita' di Jarque-Bera library(tseries) jarque.bera.test(residui) # Test di normalita' di Shapiro-Wilks shapiro.test(residui) # vediamo gli altri modelli possibili... # modello 2: Y=a+bX[,2] mod.2=lm(y~X1,data=dati) summary(mod.2) # modello 3: Y=a+bX[,3] mod.3=lm(y~X2,data=dati) summary(mod.3) # modello 0: sola intercetta mod.0=lm(y~1,data=dati) summary(mod.0) mean(y) plot(y) abline(mean(y),0,lty=3,col="red") sqrt(var(y)) # vediamo come gli R2corretti ottenuti siano inferiori al modello iniziale #---------------------------------------------------------- # * selezione delle variabili #---------------------------------------------------------- library(MASS) # Selezione bacward stepAIC(mod,direction="backward") # Selezione forward stepAIC(mod.0,scope=list(lower=mod.0,upper=mod), direction="forward",data=dati) # Selezione stepwise stepAIC(mod.0,scope=list(lower=mod.0,upper=mod),data=dati) #---------------------------------------------------------- # * multicollinearita' #---------------------------------------------------------- # Metodi per individuare la presenza dimulticollinearita' # anzitutto vediamo se il rango e' pieno rankMatrix(X)[[1]] # guardiamo alla matrice di correlazione e relativi test corr.test(dati) # calcoliamo il determinante di X'X det(t(X)%*%X) # Calcoliamo il VIF = Variance Inflation Factor library(HH) vif(dati[,-3]) #---------------------------------------------------------- # * regressione multipla + variabili dummy #---------------------------------------------------------- n=100 p=3 set.seed(123) DD=matrix(0,nrow=n,ncol=3) DD[,1]=1 DD[,2]=round(runif(n,50,100),1) DD[,3]=c(rep(0,n/2),rep(1,n/2)) b1=10 b2=4.5 b3=175 b=c(b1,b2,b3) sigma2=150 yy=round(rnorm(n,DD%*%b,sqrt(sigma2)),1) dati=data.frame(DD[,-1],yy) head(dati) names(dati)=c("X1","D","yy") boxplot(yy~factor(dati$D)) plot(density(yy[dati$D==0]),col="blue") lines(density(yy[dati$D==1]),col="pink") plot(dati$X1,yy) plot(dati$X1,yy,pch=dati$D+1) mod.dum.0=lm(yy~1,data=dati) summary(mod.dum.0) mod.dum.tutte=lm(yy~.,data=dati) summary(mod.dum.tutte) stepAIC(mod.dum.0,scope=list(lower=mod.dum.0,upper=mod.dum.tutte)) par(mfrow=c(2,2)) plot(mod.dum.tutte) #test omoschedasticit?†-eteroschedasticit?† ncvTest(mod.dum.tutte) t.test(resid(mod.dum.tutte)) shapiro.test(resid(mod.dum.tutte)) # disegnamo le due rette di regressione stimate plot(dati$X1,yy,pch=dati$D+1) abline((mod.dum.tutte$coefficients[1]), (mod.dum.tutte$coefficients[2]),lty=2,lwd=2,col="blue") abline((mod.dum.tutte$coefficients[1]+mod.dum.tutte$coefficients[3]), (mod.dum.tutte$coefficients[2]),lty=4,lwd=2,col="green") # vediamo come sarebbe il modello con la sola variabile dummy (ANOVA) mod.solo.dummy=lm(yy~D,data=dati) summary(mod.solo.dummy) mean(yy[dati$D==0]) mean(yy[dati$D==1]) mean(yy[dati$D==1])-mean(yy[dati$D==0]) mod.reg=lm(yy~dati$X1,data=dati) summary(mod.reg) plot(dati$X1,yy,pch=dati$D+1) abline(mod.reg) # Brutta cosa: osservo solo residui grandi!!! hist(resid(mod.reg)) #---------------------------------------------------------- # * dati reali (reddito) #---------------------------------------------------------- dati=read.csv("reddito.csv",header=T,sep=",") str(dati) # perche ho i factor al posto dei numeric?? # ...vai a vedere come son fatti i dati originali... dati=read.csv("reddito.csv",header=T,sep=",",dec=",") attach(dati) corr.test(dati[,-1]) corrplot(cor(dati[,-1]),method="circle") pairs(dati[,-1]) mm=lm(REDD~.,data=dati[,-1]) summary(mm) # selezione del modello library(MASS) mm0=lm(REDD~1,data=dati[,-1]) stepAIC(mm,direction="backward") stepAIC(mm0,scope=list(upper=mm,lower=mm0),direction="forward") stepAIC(mm0,scope=list(upper=mm,lower=mm0)) # scelto il modello con tutte le variabili... t.test(resid(mm)) par(mfrow=c(2,2)) plot(mm) #test omoschedasticita'-eteroschedasticita' ncvTest(mm) acf(resid(mm)) lags=1:20 sapply(lags,FUN=Box.test,x=resid(mm),type="Ljung-Box") #Comando per testare correlazione con lag da 1 a 20 shapiro.test(resid(mm)) # costruzione di una variabile dummy # partiamo dalla variabile EDU mediana=median(EDU) # dicotomizziamo la variabile EDU DUMMY.EDU=factor(EDU>mediana) boxplot(REDD~DUMMY.EDU) mm1=lm(REDD~ETA+ORE+DUMMY.EDU) summary(mm1) detach(dati) #---------------------------------------------------------- # * dati reali (circonferenza cranica) #---------------------------------------------------------- dati=read.csv("Circonferenza-Cranica.csv",header=T,sep=";") n=nrow(dati) names(dati) str(dati) # headcirc: circonferenza cranica # lenght: lunghezza / altezza del bambino # gestage: eta' gestazionale # birthwt: peso alla nascita # momage: eta' della madre al momento del parto # toxemia: presenza/assenza di tossemia della madre durante la gravidanza # rinominare le variabili con nomi piu semplici names(dati)=c("cranio","altezza","gest","peso","eta.mamma","tossemia") dati$tossemia=factor(dati$tossemia) # e attacchiamo il file in modo dal lavorare in modo piu' agevole attach(dati) # vediamo la struttura di correlazione nei dati corr.test(dati[,-6]) # cerchiamo di facilitare la lettura della matrice di correlazione library(corrplot) cor.mat=cor(dati[,-6]) corrplot(cor.mat,method="circle") pairs(dati[,-6]) boxplot(cranio~tossemia) plot(density(cranio[tossemia==0]),col="green") lines(density(cranio[tossemia==1]),col="red") mod.tutte=lm(cranio~altezza+peso+gest+eta.mamma+tossemia,data=dati) summary(mod.tutte) # metodo alternativo mod.tutte2=lm(cranio~.,data=dati) summary(mod.tutte2) # selezione del modello library(MASS) # selezione stepwise mod0=lm(cranio~1,data=dati) stepAIC(mod0,scope=list(upper=mod.tutte,lower=mod0)) # selezione backward stepAIC(mod.tutte,direction="backward") # selezione forward stepAIC(mod0,scope=list(upper=mod.tutte,lower=mod0),direction="forward") # il VIF si applica a variabili quantitative vif(dati[,-c(1,6)]) mod.finale=lm(cranio~gest+peso) summary(mod.finale) # analisi dei residui # Verifichiamo che la media sia nulla con un opportuno test t.test(resid(mod.finale)) # plottiamo i residui plot(fitted(mod.finale),resid(mod.finale),xlab=expression(hat(y)),ylab="residuals") abline(h=0) par(mfrow=c(2,2)) plot(mod.finale) #test omoschedasticita'-eteroschedasticita' ncvTest(mod.finale) bptest(mod.finale) # Per l'incorrelazione si puo' costruire il # correlogramma acf(resid(mod.finale)) # il test di Ljung-Box Box.test(resid(mod.finale),lag=1,type="Ljung-Box") Box.test(resid(mod.finale),lag=2,type="Ljung-Box") lags=1:20 sapply(lags,FUN=Box.test,x=resid(mod.finale),type="Ljung-Box") # e test di Durbin-Watson dwtest(mod.finale) shapiro.test(resid(mod.finale)) jarque.bera.test(resid(mod.finale)) hist(resid(mod.finale),prob=T,main="istogramma residui",xlab="residui") curve(dnorm(x,0,sd(resid(mod.finale))),add=T,col="red",lwd=2) # Previsione # Creaimo delle nuove osservazioni per le variabili esplicative # scelte nel modello finale, ad esempio 10 nuove osservazioni. # Occhio che dobbiamo creare un data frame # le cui variabili hanno lo stesso nome delle variabili indipendenti selezionate e compaiono # nello stesso ordine in cui sono elencate nel modello stimato. Nel nostro caso, quindi, la prima colonna # si deve chiamare gest e la seconda peso. # Inoltre, se vogliamo generare dei valori pseudo-casuali, prima fissiamo il seme così possiamo # replicare i risultati n.new=10 set.seed(9583) new.oss=data.frame(gest=sample(min(gest):max(gest),size=n.new,replace=T),peso=runif(n.new,min(peso),max(peso))) new.oss # Nell'esempio ho usato sample per generare dei valori interi ed runif per valori decimali # e ho preso il range osservato delle variabili come riferimento da cui estrarre. # Calcolo dei nuovi valori predetti per il valor medio di Y e relativi intervalli di confidenza # al 95% p_conf=predict(mod.finale,newdata=new.oss,interval="confidence") head(p_conf) # Calcolo dei nuovi valori predetti per il valore esatto di Y e relativi intervalli di previsione # al 95% p_prev=predict(mod.finale,newdata=new.oss,interval="prediction") head(p_prev) # Prima cosa da verificare e' che la prima colonna di p_conf e la prima colonna di p_pred siano uguali p_conf[,1] p_prev[,1] # La seconda e' verificare che i tali valori siano giusti. Richiamando la teoria sappiamo che la prima colonna # si ottiene sostituendo le nuove x nell'equazione stimata dell'iperpiano di regressione in modo da avere gli y.hat. # I valori predetti, dunque, non sono altro che il risultato del prodotto scalare # tra il vettore delle stime dei parametri di regressione e la # MATRICE DEL DISEGNO (attenzione: non il nuovo data frame!) # Pertanto me li calcolo a mano cos?¨ X.new=cbind(1,new.oss) # nota: ho dovuto aggiungere la colonna pari ad 1 per la matrice del disegno... X.new c(mod.finale$coefficients%*%t(X.new)) # Effettivamente sono uguali a le prime colonne viste prima! p_conf[,1] p_prev[,1] # Da notare che gli intervalli di confidenza (per il valor medio di Y) sono sempre # piu' stretti degli intervalli di previsione (per il valore esatto di Y), in quanto # per questi ultimi bisogna tener conto non solo dell'incertezza associata # alle stime dei parameteri di regressione, ma anche della stima della varianza dell'errore! # Ad esempio per la prima nuova unita' l'intervallo di confidenza e' p_conf[1,2:3] # mentre quello di previsione e' p_prev[1,2:3] # e dunque piu' largo!