#-----------------------------------------------------# # # # LABORATORIO DI METODI STATISTICI # # PER LA FINANZA # # PROF F. BARTOLUCCI # # # # LAB-07 # #-----------------------------------------------------# ########################################################################## # Lezione 7 ##########################################################################1 data = read.table("spesa_alimentare.txt",header=T) names(data) # redditor: reddito familiare # figli: numero figli # metri: metri quadri abitazione # sesso: sesso del capofamiglia (1=M, 0=F) # genitori: dummy per entrambi i genitori che lavorano (1=SI, 0=NO) # zona: zona abitazione (1=Urbana, 0=Extra-urbana) # spesa: spesa alimentare # AIM: analisi della relazione di dipendenza tra Y = spesa alimentare e le variabili esplicative rilevate # separo le variabili spesa = data$spesa reddito = data$reddito figli = data$figli metri = data$metri sesso = data$sesso genitori = data$genitori zona = data$zona # Modello di regressione lineare multipla: spesa = beta0 + beta1*reddito + beta2*figli + beta3*metri # ------------------------------------------------------------------------------------------------- mod1 = lm(spesa ~ reddito + figli + metri, data = data) summary(mod1) # estrarre gli elementi che ci interessano names(mod1) names(summary(mod1)) # La variabile metri non è significativamente diversa da zero... mod2 = lm(spesa ~ reddito + figli, data = data) summary(mod2) # variabili categoriche con due livelli # ---------------------------------------------------------- # inseriamo nel modello la variabile zona(1=Urbana, 0=Extra-Urbana) table(zona) mod3 = lm(spesa ~ reddito + figli + zona, data = data) summary(mod3) # Come cambiano i coefficienti di regressione quando inseriamo una nuova variabile nel modello? mod1$coef mod3$coef # Come interpretiamo i coefficienti associati alle variabili esplicative continue (covariate)? # Come interpretiamo i coefficienti associati alle variabili esplicative categoriche? # variabili categoriche con più di due livelli # ---------------------------------------------------------- # trasformiamo la variabile figli in una variabile categorica con 3 livelli. table(figli) # funzione factor() figli2 = factor(figli) levels(figli2) = c("nessuno", "uno", "più di uno", "più di uno") figli2 table(figli2) # inseriamo nel modello la variabile categorica figli.cat mod4 = lm(spesa ~ reddito + figli2 + zona, data = data) summary(mod4) # come interpretiamo i coefficienti associati alle variabili esplicative categoriche? # Test significatività globale anova(mod4) # Model selection # ---------------------- # procedura backward # -------------------------- # Stimo il modello completo (con tutte le covariate) mod5.1 = lm(spesa ~ reddito + figli + metri + sesso + genitori + zona) summary(mod5.1) # Eliminiamo dal modello la variabile meno significativa (con p-value max) -- sesso mod5.2 = lm(spesa ~ reddito + figli + metri + genitori + zona) summary(mod5.2) # Eliminiamo dal modello la variabile meno significativa (con p-value max) -- metri mod5.3 = lm(spesa ~ reddito + figli + genitori + zona) summary(mod5.3) # Tutte le covariate del modello sono significative al livello 0.05. # Il modello migliore mod5.3 # forward selection # -------------------------- # Stimiamo un modello per ciascuna covariata candidata alla selezione mod6.1 = lm(spesa ~ reddito) mod6.2 = lm(spesa ~ figli) mod6.3 = lm(spesa ~ metri) mod6.4 = lm(spesa ~ sesso) mod6.5 = lm(spesa ~ genitori) mod6.6 = lm(spesa ~ zona) summary(mod6.1)$adj.r.square summary(mod6.2)$adj.r.square summary(mod6.3)$adj.r.square summary(mod6.4)$adj.r.square summary(mod6.5)$adj.r.square summary(mod6.6)$adj.r.square # Fissiamo il modello a cui corrisponde l'R^2 più alto -- reddito # Aggiungiamo le rimanenti variabili al modello spesa ~ reddito mod7.1 = lm(spesa ~ reddito + figli) mod7.2 = lm(spesa ~ reddito + metri) mod7.3 = lm(spesa ~ reddito + sesso) mod7.4 = lm(spesa ~ reddito + genitori) mod7.5 = lm(spesa ~ reddito + zona) summary(mod6.1)$adj.r.square summary(mod7.1)$adj.r.square summary(mod7.2)$adj.r.square summary(mod7.3)$adj.r.square summary(mod7.4)$adj.r.square summary(mod7.5)$adj.r.square # Fissiamo il modello con variabile che comporta l'aumento più significativo dell'R^2 -- reddito + zona # Aggiungiamo le rimanenti variabili al modello spesa ~ reddito + zona mod8.1 = lm(spesa ~ reddito + zona + figli) mod8.2 = lm(spesa ~ reddito + zona + metri) mod8.3 = lm(spesa ~ reddito + zona + sesso) mod8.4 = lm(spesa ~ reddito + zona +genitori) summary(mod7.5)$adj.r.square summary(mod8.1)$adj.r.square summary(mod8.2)$adj.r.square summary(mod8.3)$adj.r.square summary(mod8.4)$adj.r.square # Fissiamo il modello con variabile che comporta l'aumento più significativo dell'R^2 -- reddito + zona + figli # Aggiungiamo le rimanenti variabili al modello spesa ~ reddito + zona + figli mod9.1 = lm(spesa ~ reddito + zona + figli + metri) mod9.2 = lm(spesa ~ reddito + zona + figli + sesso) mod9.3 = lm(spesa ~ reddito + zona + figli + genitori) summary(mod8.1)$adj.r.square summary(mod9.1)$adj.r.square summary(mod9.2)$adj.r.square summary(mod9.3)$adj.r.square # Fissiamo il modello con variabile che comporta l'aumento più significativo dell'R^2 -- reddito + zona + figli + genitori # Aggiungiamo le rimanenti variabili al modello spesa ~ reddito + zona + figli + genitori mod10.1 = lm(spesa ~ reddito + zona + figli + genitori + metri) mod10.2 = lm(spesa ~ reddito + zona + figli + genitori + sesso) summary(mod9.3)$adj.r.square summary(mod10.1)$adj.r.square summary(mod10.2)$adj.r.square # R^2 semba non aumentare in maniera consistente. Vediamo la significatività dei parametri summary(mod10.1) summary(mod10.2) # Le rimanenti variabili sono tutte non significative. Il modello migliore è spesa ~ reddito + zona + figli + genitori summary(mod9.3) # Regressione lineare multipla: multicollinearità # -------------------------------------------------------------------------- data = read.table("BloodePressure.txt", header=T) # Pt = Patient id # BP = bloode pressure (in mm Hg) # Age = eta' # Weight = peso # BSA = body surface area (in m^2) # Dur = duration of hypertension in years # Pulse = basal pulse (Pulse, in beats per minute) # Stress = stress index data = data[,-1] # AIM: verificare se esiste una relazione tra la pressione arteriosa e le variabili # età, peso, BSA, duration, pulse rate and/or stress level # Analizziamo le correlazioni tra coppie di variabili cor(data) # un po' meglio... round(cor(data), 3) # ma anche... plot(data) # cosa possiamo osservare? # Costruiamo un modello in cui la sola variabile Weight è inserita come covariata mod1 = lm(BP ~ Weight, data = data) summary(mod1) anova(mod1) # Costruiamo un modello in cui la sola variabile BSA è inserita come covariata mod2 = lm(BP ~ BSA, data = data) summary(mod2) anova(mod2) # Costruiamo un modello in cui le variabili Weight e BSA sono inserite come covariate (in questo ordine) mod3 = lm(BP ~ Weight + BSA, data = data) summary(mod3) anova(mod3) # Costruiamo un modello in cui le variabili BSA e Weight sono inserite come covariate (in questo ordine) mod4 = lm(BP ~ BSA + Weight, data = data) summary(mod4) anova(mod4) # 1° effetto: stime non robuste dei parametri coef(mod1) coef(mod2) coef(mod3) # Cambiando le variabili inserite nel modello si ottengono risultati molto differenti! # 2° effetto: aumento degli se round(summary(mod1)$coeff, 3) round(summary(mod2)$coeff, 3) round(summary(mod4)$coeff, 3) # cambiano le conclusioni sulla significatività dei parametri! # 3° effetto: la quota di devianza spiegata da ciascun predittore cambia a # seconda dei predittori già presenti nel modello anova(mod1)[2] anova(mod2)[2] anova(mod3)[2] anova(mod4)[2] # verificare la presenza di multicollinearità # ------------------------------------------------ # Hp: siamo interessati al modello che mette in relazione BP con Weight, BSA Pulse e Age mod = lm(BP ~ Weight + BSA + Pulse + Age, data = data) summary(mod) # VIF - Variance inflation factor # Calcolo VIF1 # Regredisco Weight su BSA Pulse e Age mod1 = lm(Weight ~ BSA + Pulse + Age, data = data) summary(mod1) # R^2 aggiustato: 1- ((RSS/n-p-1)/(RSS/n-1)) r1 = summary(mod1)$adj.r.squared VIF1 = 1/(1-r1) # Calcolo VIF2 # Regredisco BSA su Weighe Pulse e Age mod2 = lm(BSA ~ Weight + Pulse + Age, data = data) summary(mod2) # R^2 aggiustato: 1- ((RSS/n-p-1)/(RSS/n-1)) r2 = summary(mod2)$adj.r.squared VIF2 = 1/(1-r2) # Calcolo VIF3 # Regredisco Pulse su Weight, BSA e Age mod3 = lm(Pulse ~ Weight + BSA + Age, data = data) summary(mod3) # R^2 aggiustato: 1- ((RSS/n-p-1)/(RSS/n-1)) r3 = summary(mod3)$adj.r.squared VIF3 = 1/(1-r3) # Calcolo VIF4 # Regredisco Age su Weight, BSA e Pulse mod4 = lm(Age ~ Weight + BSA + Pulse, data = data) summary(mod4) # R^2 aggiustato: 1- ((RSS/n-p-1)/(RSS/n-1)) r4 = summary(mod4)$adj.r.squared VIF4 = 1/(1-r4) VIFk = matrix(c(VIF1, VIF2, VIF3, VIF4), 1,4) colnames(VIFk) = c("Weight", "BSA", "Pulse", "Age") VIFk # Lo se dei coefficienti Weight e BSA è inflazionata di oltre 4 punti a causa della correlazione con almeno # uno dei rimanenti predittori # VIF medio VIFm = (VIF1 + VIF2 + VIF3 + VIF4)/4 VIFm # più velocemente.... library(car) vif(mod) # proviamo ad eliminare la variabile Weight dal modello mod = lm(BP ~ BSA + Pulse + Age, data = data) vif(mod) summary(mod) # Regressione multivariata # ------------------------- data = read.table("regr_multivariata.txt",header=T,sep="\t") str(data) colnames(data) = c("eta", "sesso", "titolo", "reddito", "anni") # stimiamo il modello per il reddito mod1 = lm(reddito ~ eta + sesso + titolo, data = data) summary(mod1) # stimiamo il modello per gli anni di lavoro mod2 = lm(anni ~ eta + sesso + titolo, data) summary(mod2) # estraiamo la matrice del disegno X = model.matrix(mod1) # oppure: X = model.matrix(mod2) head(X) # creiamo la matrice delle variabili risposta Y = cbind(data$reddito,data$anni) head(Y) colnames(Y) = c("reddito", "anni") # calcoliamo i coefficienti di regressione Bh = solve(t(X)%*%X)%*%t(X)%*%Y Bh # calcoliamo i valori previsti delle due variabili risposta Yh = X%*%Bh head(Yh) # Calcoliamo i residui Eh = Y-Yh head(Eh) # Calcoliamo la matrice di varianze e covarianze tra le variabili risposta n = nrow(Yh) npar = dim(X)[2] Sigmah = t(Eh) %*% Eh/ (n-npar) Sigmah # Correlazione tra reddito ed anni di lavoro Cov = Sigmah[1,2] Cov Sigma1 = sqrt(Sigmah[1,1]) Sigma1 Sigma2 = sqrt(Sigmah[2,2]) Sigma2 Corr = Cov/(Sigma1*Sigma2) Corr # Regressione logistica # ------------------------------ # lettura dati birthwt = read.table("birthwt.txt", header=TRUE, quote="\"", stringsAsFactors=FALSE) # descrizione del ds # LOW = Low Birth Weight (0 = Birth Weight >= 2500g, 1 = Birth Weight < 2500g) # AGE = Age of the Mother in Years # KG = Weight of the mother in KG at the Last Menstrual Period # RACE = Race of the mother (1 = White, 2 = Black, 3 = Other) # FVT = Number of Physician Visits During the First Trimester (0 = None, 1 = One, 2 = Two, etc.) ## trasformiamo race in factor birthwt$race = factor(birthwt$race) levels(birthwt$race) = c("White","Black","Other") # partiamo dal modello con la sola variabile kg mod1 = glm(cbind(low, 1-low)~kg , family = "binomial", data = birthwt) summary(mod1) # il rischio di evento su scala logaritmica diminuisce di 0.03 quando aumenta di una unità il peso della madre ## costruiamo l'odds ratio per interpretare beta1 exp(coef(mod1)[2]) # Aumentando di una unità il peso della madre l'odds ratio è pari a 0.97 # inseriamo la variabile race (fattore) mod2 = glm(cbind(low, 1-low)~kg + factor(race), family = "binomial", data = birthwt) summary(mod2) # commenti? # l'intercetta non è significativamente diversa da 0. # Il successo (bambino sottopeso) e l'insuccesso (bambino normopeso) sono egualmente probabili # Aumentando il peso della madre di una unità (un KG), la probabilità di avere un figlio sottopeso diminuisce (segno negativo) # Quando aumenta di una unità (un kg) il peso della madre al concepimento, # il rischio di evento (sottopeso) diminuisce 0.03 su scala logaritmica # La categoria "Withe" per la variabile RACE è utilizzata come categoria di riferimento # Passando da "White" a "Black" si osserva un aumento della probabilità di successo di 1.08 su scala logaritmica. # Non sembra esserci una differenza significativa tra "White" e "Others" # odds ratio per KG e RACE OR = exp(coef(mod3)[-c(1, 4)]) OR ## Quando aumenta di una unità il peso della madre al concepimento, il rischio di evento diminuisce di 3 punti percentuali ## Il rischio di avere un bambino sottopeso è 2.95 più elevato per i neri rispetto ai bianchi e alle altre razze # Stimo gli intervalli di confidenza dei coefficienti di regressione confint(mod2) # Per facilitare l'interpretazione del modello è utile calcolare l'esponenziale dei coefficienti di regressione, così da # interpretarli come odds-ratio: exp(confint(mod2)) # inserisco nel modello la variabile age mod3 = glm(cbind(low, 1-low)~kg + factor(race) + age, family = "binomial", data = birthwt) summary(mod3) # il coefficiente associato a age non è statisticamente significativo (Test di Wald). # Confrontiamo con il LRT i modelli con e senza la variabile age anova(mod2, mod3, test = "Chisq") # Il test suggerisce (come atteso) una non significatività del coefficiente Age. # Calcolo le stime delle probabilità individuali di # essere ammessi ad una graduate school prob = predict(mod2, type= "response") prob