%%%%%% TENTATIVI DI REGRESSIONE IB <- read.table(file="indicatori_benessere.txt",header=TRUE) biplot(princomp(IB)) PLIC<-IB$PLIC SC<-IB$SC SA<-IB$SA.SC TD<-IB$TD TMI<-IB$TMI % regressione semplice: plot(TD,TMI) plot(IB) cor(IB) % produce: > cor(IB) PLIC SC SA.SC TD TMI PLIC 1.0000000 0.3224224 -0.4109997 -0.3665858 -0.4435252 SC 0.3224224 1.0000000 -0.8418786 -0.8502698 -0.4835062 SA.SC -0.4109997 -0.8418786 1.0000000 0.9054330 0.5138328 TD -0.3665858 -0.8502698 0.9054330 1.0000000 0.4869491 TMI -0.4435252 -0.4835062 0.5138328 0.4869491 1.0000000 SA.TD <- lm(SA ~ TD) SA.TD % produce: > SA.TD Call: lm(formula = SA ~ TD) Coefficients: (Intercept) TD -1.405e-10 9.054e-01 plot(TD,SA) x <- ((min(TD)*10):(max(TD)*10))/10 y <- SA.TD$coefficient[1]+SA.TD$coefficient[2]*x lines(x,y) % plot alternativo: plot(TD,SA) abline(SA.TD$coefficient) res<-SA.TD$residuals delta <- sd(res)*qnorm(1-0.05/2) abline(SA.TD$coefficient+c(delta,0), lty=2) abline(SA.TD$coefficient+c(-delta,0), lty=2) $ giudichiamo la bontà della regressione: summary(SA.TD) % produce: > summary(SA.TD) Call: lm(formula = SA ~ TD) Residuals: Min 1Q Median 3Q Max -0.5344 -0.3521 -0.1114 0.3433 0.8596 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.405e-10 9.752e-02 -1.44e-09 1 TD 9.054e-01 1.001e-01 9.05 4.06e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.4361 on 18 degrees of freedom Multiple R-Squared: 0.8198, Adjusted R-squared: 0.8098 F-statistic: 81.89 on 1 and 18 DF, p-value: 4.058e-08 cor(TD,SA)^2 % è R^2. $ NOTA: p-value piccolo significa che il test sarebbe stato significativo anche con elevata significatività, quindi si rigetta l'ipotesi nulla. Nel caso di un singolo coefficiente, l'ipotesi nulla è che il modello senza di esso produca risultati comparabili a quello con esso. Se viene rifiutata, significa che il coefficiente è importante. Se non viene rifiutata, significa che non abbiamo buone ragioni per ritenere importante quel coefficiente. % Analisi dei residui: plot(SA.TD) % confrontare il primo grafico con: plot(TD,SA) x <- ((min(TD)*10):(max(TD)*10))/10 y <- SA.TD$coefficient[1]+SA.TD$coefficient[2]*x lines(x,y) % nel seguente modo si estraggono i residui: res<-SA.TD$residuals plot(res) % ma il plot è contro l'indice. % Regressione multipla: mod1 <- lm(TMI ~ PLIC + TD) mod1 % produce: > mod1 Call: lm(formula = TMI ~ PLIC + TD) Coefficients: (Intercept) PLIC TD 5.845e-11 -3.062e-01 3.747e-01 % notare il segno dei coefficienti summary(mod1) % produce: > summary(mod1) Call: lm(formula = TMI ~ PLIC + TD) Residuals: Min 1Q Median 3Q Max -1.52928 -0.62822 0.01491 0.52827 1.45582 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.845e-11 1.952e-01 2.99e-10 1.0000 PLIC -3.062e-01 2.152e-01 -1.422 0.1730 TD 3.747e-01 2.152e-01 1.741 0.0998 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.8729 on 17 degrees of freedom Multiple R-Squared: 0.3183, Adjusted R-squared: 0.2381 F-statistic: 3.968 on 2 and 17 DF, p-value: 0.03853 % il modello è scadente plot(mod1) % interpretare i residui. %%%%% FATTORI NASCOSTI PCA <- princomp(IB) plot(PCA) summary(PCA) % da entrambi emerge dimensione (numero di variabili nascoste): 2-3. % Se vogliamo un gruppo uni-dimensionale: IB.new <- matrix(nrow=20,ncol=3) IB.new[,1]<-SC IB.new[,2]<-SA IB.new[,3]<-TD PCA.new<-princomp(IB.new) plot(PCA.new) % tentativo sbagliato: IB.new2 <- matrix(nrow=20,ncol=3) IB.new2[,1]<-PLIC IB.new2[,2]<-SA IB.new2[,3]<-TMI PCA.new2<-princomp(IB.new2) plot(PCA.new2) % Altro modo per scartare variabili e ridursi a 1 dim: PCA$loadings % produce: > PCA$loadings Loadings: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 PLIC -0.310 0.769 -0.553 SC -0.491 -0.309 -0.813 SA.SC 0.512 0.216 0.120 -0.433 -0.699 TD 0.506 0.279 0.115 -0.381 0.713 TMI 0.379 -0.435 -0.816 Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 SS loadings 1.0 1.0 1.0 1.0 1.0 Proportion Var 0.2 0.2 0.2 0.2 0.2 Cumulative Var 0.2 0.4 0.6 0.8 1.0 % e qui si scelgono le variabili più legate a comp1. % Per metterle ancor più in evidenza: varimax(PCA$loadings[,1:3]) % produce: > varimax(PCA$loadings[,1:3]) $loadings Loadings: Comp.1 Comp.2 Comp.3 PLIC 0.996 SC -0.577 SA.SC 0.566 TD 0.589 TMI -1.000 Comp.1 Comp.2 Comp.3 SS loadings 1.0 1.0 1.0 Proportion Var 0.2 0.2 0.2 Cumulative Var 0.2 0.4 0.6 $rotmat [,1] [,2] [,3] [1,] 0.8721064 -0.3097967 -0.3787566 [2,] 0.4632945 0.7718482 0.4354408 [3,] 0.1574444 -0.5552265 0.8166608 % varimax(PCA$loadings) non serve a niente. % Solo due comp non separano: varimax(PCA$loadings[,1:2]) %%%% Che modello lineare sta dietro questi numeri???