## a.a. 24 25 ##code per Regressione base ##1 Ottobre 2024 ## Statistica Aziendale setwd("/Users/mariagraziapittau/Library/CloudStorage/Dropbox-Universita'DegliStudiDiRomaLaSapienza/grazia.pittau@uniroma1.it/Lezioni/Statistica_Aziendale/R/Regressione") ## https://hbr.org/2015/11/a-refresher-on-regression-analysis library(foreign) library(arm) kidiq <- read.dta(file="kidiq.dta") kidiq[1:5,] ### 1.Stima e sintesi del modello di Regressione in R--un SOLO predittore BINARIO fit.1 <- lm (kid_score ~ mom_hs, data=kidiq) display(fit.1) summary(fit.1) ##il coeff. nel caso di regressione con input binario ##e' pari alla differenza tra le medie delle due sottopopolazioni Media1<-mean(kidiq[, "kid_score"][kidiq[,"mom_hs"]==1]) Media0<-mean(kidiq[, "kid_score"][kidiq[,"mom_hs"]==0]) Media1-Media0 ### Rappresentazione grafica dei dati e del modello stimato plot (jitter(kidiq$mom_hs, 0.3), jitter(kidiq$kid_score, 0.3), xlab="Madri con diploma superiore", ylab="Punteggio dei figli", xaxt="n") curve (coef(fit.1)[1] + coef(fit.1)[2]*x, add=TRUE) yrs <- c(1, 2) axis(1, yrs, tick = T, col.axis = 'black') ### 2.Stima e sintesi del modello di Regressione in R--un SOLO predittore CONTINUO fit.2 <- lm (kid_score ~ mom_iq, data=kidiq) display(fit.2) summary(fit.2) ### Rappresentazione grafica dei dati e del modello stimato plot (kidiq$mom_iq, kidiq$kid_score, xlab="Punteggio IQ delle madri", ylab="Punteggio dei figli") curve (coef(fit.2)[1] + coef(fit.2)[2]*x, add=TRUE) ### 3.Stima e sintesi del modello di Regressione in R-- predittori MULTIPLI fit.3 <- lm (kid_score ~ mom_hs + mom_iq, data=kidiq) display(fit.3) print(fit.3) summary(fit.3) ### Rappresentazione grafica dei dati e del modello stimato ## oss. abbiamo due rette parallele con stessa pendenza colors <- ifelse (kidiq$mom_hs==1, "red", "blue") plot (kidiq$mom_iq, kidiq$kid_score, xlab="Punteggio IQ delle madri", ylab="Punteggio dei figli", col=colors, pch=20) curve (cbind (1, 1, x) %*% coef(fit.3), add=TRUE, col="red") curve (cbind (1, 0, x) %*% coef(fit.3), add=TRUE, col="blue") # Una sequenza alternativa di comandi in R plot(kidiq$mom_iq, kidiq$kid_score,xlab="Punteggio IQ delle madri",ylab="Punteggio dei figli",type="n") points (kidiq$mom_iq[kidiq$mom_hs==1], kidiq$kid_score[kidiq$mom_hs==1], pch=20, col="gray") points (kidiq$mom_iq[kidiq$mom_hs==0], kidiq$kid_score[kidiq$mom_hs==0], pch=20, col="black") curve (cbind (1, 1, x) %*% coef(fit.3), add=TRUE, col="gray") curve (cbind (1, 0, x) %*% coef(fit.3), add=TRUE, col="black") ##4.Stima e sintesi del modello di Regressione in R-- MODELLO CON INTERAZIONE fit.4a <- lm (kid_score ~ mom_hs + mom_iq + mom_hs*mom_iq, data=kidiq) display(fit.4a) ### Rappresentazione grafica dei dati e del modello stimato par (mfrow=c(1,2)) colors <- ifelse (kidiq$mom_hs==1, "red", "black") plot (kidiq$mom_iq, kidiq$kid_score, xlab="Punteggio IQ delle madri", ylab="Punteggio dei figli", col=colors, pch=20, xlim=c(0, 140), ylim=c(-50, 140)) curve (cbind (1, 1, x, 1*x) %*% coef(fit.4a), add=TRUE, col="red") curve (cbind (1, 0, x, 0*x) %*% coef(fit.4a), add=TRUE, col="black") colors <- ifelse (kidiq$mom_hs==1, "red", "black") mom_iqb <- cbind(0, kidiq$mom_iq) kid_scoreb <- cbind(0, kidiq$kid_score) plot (range(mom_iqb), range(kid_scoreb), xlab="Punteggio IQ delle madri", ylab="Punteggio dei figli", col=colors, pch=20, type="n") points (kidiq$mom_iq, kidiq$kid_score, col=colors, pch=20) curve (cbind (1, 1, x, 1*x) %*% coef(fit.4a), add=TRUE, col="red") curve (cbind (1, 0, x, 0*x) %*% coef(fit.4a), add=TRUE, col="black") legend(20, 100, c("HS", "No_HS"), c("red","black"),cex=0.5) ## Modello con interazione fumo e radon par (mar=c(5,5,4,1)+.1) plot (c(0,12.5),c(0,.25), type="n", xaxs="i", yaxs="i", xlab="Esposizione al radon all'interno di una casa (pCi/L)", ylab="Probabilita' di contrarre tumore al polmone", main="Esempio di interazione", cex.main=1, cex.axis=1, cex.lab=1) lines (c(0,20),.07409+c(0,20)*.0134) lines (c(0,20),.00579+c(0,20)*.0026) text (10, .07409+10*.0134 - .02, "Fumatori", cex=1) text (10, .00579+10*.0026 + .01, "Non-fumatori", cex=1) ##Variabili categoriche ##https://hbswk.hbs.edu/item/kids-of-working-moms-grow-into-happy-adu lts ##mom.work = 1: se la madre non ha lavorato nei primi tre anni di vita del bambino ## mom.work = 2: se la madre ha lavorato nel secondo o terzo anno di vita del bambino ##mom.work = 3: se la madre ha lavorato part-time nel primo anno di vita del bambino ## mom.work = 4: se la madre ha lavorato a tempo pieno nel primo anno di vita del bambini kidiq$mom_work <- as.factor(kidiq$mom_work) fit.5 <- lm (kid_score ~ mom_iq+mom_work+mom_hs, data=kidiq) display(fit.5) kidiq$mom_work <- relevel(kidiq$mom_work,"4") fit.5a <- lm (kid_score ~mom_hs+mom_iq+mom_work+mom_iq:mom_work, data=kidiq) display(fit.5a) ##Analisi dei residui ## Considero il modello fit.2 fit.2 <- lm (kid_score ~ mom_iq, data=kidiq) residui <- resid(fit.2) residui <- fit.2$residual y.hat<-fit.2$fitted.values library(ggplot2) ## Valori previsti versus valori osservati qplot (kidiq$kid_score, y.hat, xlim=c(40, 150), ylim=c(40, 150), xlab="Punteggi osservati", ylab="Punteggi stimati", main="Valori stimati verso valori previsti", alpha=I(1/5)) + geom_abline(intercept=0, slope=1, linetype="dotted") + theme_bw() ## Valori previsti e residui qplot (y.hat, residui, xlab="Punteggi stimati", ylab="Residui", main="Residui verso valori previsti", alpha=I(1/5)) + geom_abline(intercept=0, slope=0, linetype="dotted") + theme_bw() ## Valori della predittore e residui plot (kidiq$mom_iq, residui, xlab="Punteggio IQ delle madri", ylab="residui") abline(h=0) abline(h=-sigma.hat(fit.2), lty=2) abline(h=+sigma.hat(fit.2), lty=2) ##Istogramma dei residui hist(residui, breaks=20, freq=F, ylab="densita'", main="Istogramma dei residui") ##su frequenze relative, altrimenti ho il problema della scala lines(density(residui), col="red") ### Prediction x.new <- data.frame (mom_hs=1, mom_iq=100) predict (fit.3, x.new, interval="prediction", level=0.95) ##Plot dei coefficienti coefplot(fit.5, var.las=1, mar=c(1,8,5.1,2)) mtext( bquote( ( R^2==.(round(display(fit.5)$r.squared, 2)) ) ) ) ############################### ###KERNEL ESTIMATION ############################### residui <- fit.2$residual hist(residui, nclass=60) par(mfrow=c(2,1)) ##kernel ## kernel densities with different kernel functions plot(density(residui, n=length(residui), bw = bw.nrd(residui), kernel="gaussian"), xlab="Residui", ylab="Density", main="Kernel densities with different kernel functions") lines(density(residui, n=length(residui), bw = bw.nrd(residui), kernel="epanechnikov"), col="red") lines(density(residui, n=length(residui), bw = bw.nrd(residui), kernel="triangular"), col="blue") lines(density(residui, n=length(residui), bw = bw.nrd(residui), kernel="biweight"), col="green") legend(25, 0.1, c("Gaussian", "Epanechnikov", "Triangular", "Biweight"), col=c("black", "red","blue", "green"), pch=19, cex=0.8) ##Different bandwidths ##nrd Normal Reference Distribution plot(density(residui, n=length(residui), bw = bw.SJ(residui), kernel="gaussian"), xlab="Residui", ylab="Density", main="Kernel densities with different bandwidths") lines(density(residui, n=length(residui), bw = bw.nrd0(residui), kernel="gaussian"), col="red") lines(density(residui, n=length(residui), bw = bw.nrd(residui), kernel="gaussian"), col="blue") legend(25, 0.1, c("Sheather & Jones", "Silverman's of rule of thumb", "Scott with factor 1.06"), col=c("black", "red","blue"), pch=19, cex=0.7)