Prime analisi

Analizzeremo dei dati medici riguardanti il cancro alla prostata, per prima cosa caricheremi i packages di cui avremo bisogno, i dati e inizieremo con le prima analisi.

library(ggplot2)
library(corrplot)
## corrplot 0.92 loaded
data<-read.csv("prostate.csv",header=TRUE,row.names = 1)
data<-data[,-10]
summary(data)
##      lcavol           lweight           age             lbph        
##  Min.   :-1.3471   Min.   :2.375   Min.   :41.00   Min.   :-1.3863  
##  1st Qu.: 0.5128   1st Qu.:3.376   1st Qu.:60.00   1st Qu.:-1.3863  
##  Median : 1.4469   Median :3.623   Median :65.00   Median : 0.3001  
##  Mean   : 1.3500   Mean   :3.629   Mean   :63.87   Mean   : 0.1004  
##  3rd Qu.: 2.1270   3rd Qu.:3.876   3rd Qu.:68.00   3rd Qu.: 1.5581  
##  Max.   : 3.8210   Max.   :4.780   Max.   :79.00   Max.   : 2.3263  
##       svi              lcp             gleason          pgg45       
##  Min.   :0.0000   Min.   :-1.3863   Min.   :6.000   Min.   :  0.00  
##  1st Qu.:0.0000   1st Qu.:-1.3863   1st Qu.:6.000   1st Qu.:  0.00  
##  Median :0.0000   Median :-0.7985   Median :7.000   Median : 15.00  
##  Mean   :0.2165   Mean   :-0.1794   Mean   :6.753   Mean   : 24.38  
##  3rd Qu.:0.0000   3rd Qu.: 1.1787   3rd Qu.:7.000   3rd Qu.: 40.00  
##  Max.   :1.0000   Max.   : 2.9042   Max.   :9.000   Max.   :100.00  
##       lpsa        
##  Min.   :-0.4308  
##  1st Qu.: 1.7317  
##  Median : 2.5915  
##  Mean   : 2.4784  
##  3rd Qu.: 3.0564  
##  Max.   : 5.5829
corrplot(cor(data),method= "number")

Da una prima seplice analisi possiamo subito notare come le correlazioni siano tutte positive rispetto alla nostra variabile risposta lpsa. La variabile indipendente che sembra avere la maggior correlazione con lpsa è lcavol, mentre age e lbph sono molto vicine all’incorrelazione.

Dal summary() possiamo vedere che si tratta di tutte variabili numeriche e distribuite in maniera uniforme, l’unica che potrebbe destare qualche dubbio è pgg45 ma applicare trasformazioni a dati medici non è sempre banale, quindi la terremo così.

Polinomiale vs Lineare

Stimiamo un modello lineare e poi uno o più polinomiali per metterli a confronto.

Inizialmente faremo modelli contnenti tutte le var.indipendenti, poi man mano proveremo a creare modelli più compatti per valutarne l’effetto sul confronto.

m<-lm(lpsa~.,data)

mpol2<-lm(lpsa~ poly(data$lcavol,2,raw=T)+poly(data$lweight,2,raw=T)+
            poly(data$age,2,raw=T)+poly(data$lbph,2,raw=T)+poly(data$svi,2,raw=T)+poly(data$lcp,2,raw=T)+
            poly(data$gleason,2,raw=T)+poly(data$pgg45,2,raw=T),data )

mpol3<-lm(lpsa~ poly(data$lcavol,3,raw=T)+poly(data$lweight,3,raw=T)+
            poly(data$age,3,raw=T)+poly(data$lbph,3,raw=T)+poly(data$svi,3,raw=T)+poly(data$lcp,3,raw=T)+
            poly(data$gleason,3,raw=T)+poly(data$pgg45,3,raw=T),data )

Essendo modelli nidificati possiamo confrontarli con la funzione anova()

anova(m,mpol2,mpol3)

Da un primo confronto sembra che i 2 modelli polinomiali non offrano un vantaggio tale da giustificare la loro implementazione, il modello lineare sembra più che sufficiente per i dati a disposizione.

summary(m)
## 
## Call:
## lm(formula = lpsa ~ ., data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.76644 -0.35510 -0.00328  0.38087  1.55770 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.181561   1.320568   0.137  0.89096    
## lcavol       0.564341   0.087833   6.425 6.55e-09 ***
## lweight      0.622020   0.200897   3.096  0.00263 ** 
## age         -0.021248   0.011084  -1.917  0.05848 .  
## lbph         0.096713   0.057913   1.670  0.09848 .  
## svi          0.761673   0.241176   3.158  0.00218 ** 
## lcp         -0.106051   0.089868  -1.180  0.24115    
## gleason      0.049228   0.155341   0.317  0.75207    
## pgg45        0.004458   0.004365   1.021  0.31000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6995 on 88 degrees of freedom
## Multiple R-squared:  0.6634, Adjusted R-squared:  0.6328 
## F-statistic: 21.68 on 8 and 88 DF,  p-value: < 2.2e-16
summary(mpol2)
## 
## Call:
## lm(formula = lpsa ~ poly(data$lcavol, 2, raw = T) + poly(data$lweight, 
##     2, raw = T) + poly(data$age, 2, raw = T) + poly(data$lbph, 
##     2, raw = T) + poly(data$svi, 2, raw = T) + poly(data$lcp, 
##     2, raw = T) + poly(data$gleason, 2, raw = T) + poly(data$pgg45, 
##     2, raw = T), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.90917 -0.35545 -0.01612  0.36125  1.72973 
## 
## Coefficients: (1 not defined because of singularities)
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      2.0134874  9.8766226   0.204 0.838971    
## poly(data$lcavol, 2, raw = T)1   0.5193401  0.1443082   3.599 0.000549 ***
## poly(data$lcavol, 2, raw = T)2   0.0229237  0.0559016   0.410 0.682836    
## poly(data$lweight, 2, raw = T)1 -1.4177739  2.0515076  -0.691 0.491486    
## poly(data$lweight, 2, raw = T)2  0.3015572  0.2829369   1.066 0.289676    
## poly(data$age, 2, raw = T)1     -0.1210170  0.1087979  -1.112 0.269295    
## poly(data$age, 2, raw = T)2      0.0007981  0.0008929   0.894 0.374048    
## poly(data$lbph, 2, raw = T)1     0.1345019  0.0640065   2.101 0.038716 *  
## poly(data$lbph, 2, raw = T)2    -0.1488016  0.0827317  -1.799 0.075806 .  
## poly(data$svi, 2, raw = T)1      0.7252933  0.2553029   2.841 0.005688 ** 
## poly(data$svi, 2, raw = T)2             NA         NA      NA       NA    
## poly(data$lcp, 2, raw = T)1     -0.1175504  0.0961782  -1.222 0.225171    
## poly(data$lcp, 2, raw = T)2     -0.0434383  0.0615019  -0.706 0.482033    
## poly(data$gleason, 2, raw = T)1  1.4493619  2.1593815   0.671 0.504007    
## poly(data$gleason, 2, raw = T)2 -0.0966872  0.1465923  -0.660 0.511404    
## poly(data$pgg45, 2, raw = T)1    0.0063049  0.0153963   0.410 0.683248    
## poly(data$pgg45, 2, raw = T)2   -0.0000273  0.0001813  -0.151 0.880661    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.699 on 81 degrees of freedom
## Multiple R-squared:  0.6906, Adjusted R-squared:  0.6334 
## F-statistic: 12.06 on 15 and 81 DF,  p-value: 5.942e-15
summary(mpol3)
## 
## Call:
## lm(formula = lpsa ~ poly(data$lcavol, 3, raw = T) + poly(data$lweight, 
##     3, raw = T) + poly(data$age, 3, raw = T) + poly(data$lbph, 
##     3, raw = T) + poly(data$svi, 3, raw = T) + poly(data$lcp, 
##     3, raw = T) + poly(data$gleason, 3, raw = T) + poly(data$pgg45, 
##     3, raw = T), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.75410 -0.39025 -0.02384  0.22693  1.72016 
## 
## Coefficients: (2 not defined because of singularities)
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      5.889e+01  1.719e+02   0.343 0.732792    
## poly(data$lcavol, 3, raw = T)1   5.339e-01  1.468e-01   3.636 0.000508 ***
## poly(data$lcavol, 3, raw = T)2  -9.348e-02  1.410e-01  -0.663 0.509324    
## poly(data$lcavol, 3, raw = T)3   3.526e-02  3.810e-02   0.926 0.357649    
## poly(data$lweight, 3, raw = T)1 -1.598e+01  1.533e+01  -1.042 0.300752    
## poly(data$lweight, 3, raw = T)2  4.407e+00  4.280e+00   1.030 0.306468    
## poly(data$lweight, 3, raw = T)3 -3.801e-01  3.933e-01  -0.967 0.336870    
## poly(data$age, 3, raw = T)1     -7.745e-02  1.030e+00  -0.075 0.940264    
## poly(data$age, 3, raw = T)2     -1.907e-04  1.724e-02  -0.011 0.991207    
## poly(data$age, 3, raw = T)3      7.184e-06  9.490e-05   0.076 0.939857    
## poly(data$lbph, 3, raw = T)1     2.766e-01  2.337e-01   1.183 0.240428    
## poly(data$lbph, 3, raw = T)2    -7.202e-02  1.481e-01  -0.486 0.628296    
## poly(data$lbph, 3, raw = T)3    -6.176e-02  1.012e-01  -0.610 0.543629    
## poly(data$svi, 3, raw = T)1      7.913e-01  2.657e-01   2.978 0.003919 ** 
## poly(data$svi, 3, raw = T)2             NA         NA      NA       NA    
## poly(data$svi, 3, raw = T)3             NA         NA      NA       NA    
## poly(data$lcp, 3, raw = T)1     -3.227e-01  1.864e-01  -1.731 0.087598 .  
## poly(data$lcp, 3, raw = T)2     -1.851e-01  1.054e-01  -1.756 0.083225 .  
## poly(data$lcp, 3, raw = T)3      8.314e-02  5.581e-02   1.490 0.140531    
## poly(data$gleason, 3, raw = T)1 -1.603e+01  7.195e+01  -0.223 0.824318    
## poly(data$gleason, 3, raw = T)2  2.372e+00  9.927e+00   0.239 0.811828    
## poly(data$gleason, 3, raw = T)3 -1.130e-01  4.505e-01  -0.251 0.802685    
## poly(data$pgg45, 3, raw = T)1   -3.721e-02  3.637e-02  -1.023 0.309650    
## poly(data$pgg45, 3, raw = T)2    1.156e-03  8.892e-04   1.301 0.197448    
## poly(data$pgg45, 3, raw = T)3   -8.706e-06  6.406e-06  -1.359 0.178284    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7006 on 74 degrees of freedom
## Multiple R-squared:  0.716,  Adjusted R-squared:  0.6316 
## F-statistic: 8.481 on 22 and 74 DF,  p-value: 1.209e-12

Attraverso i summary() vediamo i coefficienti più significativi per ogni modello e proviamo a crearne altri 3 che contengano solo questi.

m2<-lm(lpsa~ lcavol+lweight+svi,data)
mpol2.2<-lm(lpsa~ I(lcavol^2)+ lbph+ svi,data)
mpol3.2<-lm(lpsa~ I(lcavol^2)+ svi,data)

BIC(m2,mpol2.2,mpol3.2)

Confrontando i 3 modelli tramite il calcolo del BIC notiamo che, anche dopo aver scelto i modelli più performanti, quello lineare sembra presentare qualità migliori.

Stima e visualizzazione di un modello polinomiale

Per fare ciò useremo ggplot2 che consente in poche righe di codice di fare entrambe le cose contemporaneamente.

ggplot(data,aes(x=lcavol,y=lpsa))+
  geom_point(alpha=0.7)+
  geom_smooth(method="lm",formula=y ~ x,se=F)+
  geom_smooth(method="lm", formula = y ~ poly(x,2,raw=T), col="red",se=F)+
  geom_smooth(method="lm", formula = y ~ poly(x,3,raw=T), col="green",se=F)+
  geom_smooth(method="lm", formula = y ~ poly(x,4,raw=T), col="purple",se=F)+
  geom_smooth(method="lm", formula = y ~ poly(x,10,raw=T), col="orange",se=F)

Forest Plot

Un modo per visualizzare in manier semplice ed immediata gli effetti delle covariate sulla variabile risposta è quello di usare un forest plot: in questo grafico vengono rappresentati i coefficienti di ogni variabile indipendente con i relativi standard error; più la variabile è significativa nel modello e più il suo coefficiente all’interno del forest plot sarà lontano dallo 0.

library(arm) 
## Caricamento del pacchetto richiesto: MASS
## Caricamento del pacchetto richiesto: Matrix
## Caricamento del pacchetto richiesto: lme4
## 
## arm (Version 1.13-1, built: 2022-8-25)
## Working directory is C:/Uni_magistrale/Apprendimento_statistico/R/Laboratorio_3
## 
## Caricamento pacchetto: 'arm'
## Il seguente oggetto è mascherato da 'package:corrplot':
## 
##     corrplot
par(mar=c(5,7,4,4))
coefplot(m)                

coefplot(mpol2)

In questo caso abbianmo creato i forest plot per il modello lineare e per il omdello polinomiale di grado 2, non ci sono molte sorprese per quanto riguarda le variabili significative e la loro significatività, infatti abbiamo le stesse informazioni che abbiamo potuto osservare prima attraverso i summary() dei 2 modelli, ma la visualizzazione e l’interpretazione è sicuramente più veloce e chiara.

Step function

Stimiamo delle regressioni a gradini che differiscano per numero di break points e per variabile indipemndente utilizzata.

Nel primo caso useremo age come variabili indipendente:

cut<-cut(sort(data$age),4)
mcut<-lm(lpsa~cut,data)
pred<-predict(mcut,data)

cut2<-cut(sort(data$age),10)
mcut2<-lm(lpsa~cut2,data)
pred2<-predict(mcut2,data)

cut3<-cut(sort(data$age),2)
mcut3<-lm(lpsa~cut3,data)
pred3<-predict(mcut3,data)

cut4<-cut(sort(data$age),100)
mcut4<-lm(lpsa~cut4,data)
pred4<-predict(mcut4,data)

par(mfrow=c(2,2))
plot(data$age,data$lpsa,cex=0.7,main="4 Break points",xlab="Age",ylab="lpsa")
 lines(sort(data$age),pred,type="p",col="red",pch=15,cex=0.7)
plot(data$age,data$lpsa,cex=0.7,main="10 Break points",xlab="Age",ylab="lpsa") 
 lines(sort(data$age),pred2,type="p",col="green",pch=15,cex=0.7)
plot(data$age,data$lpsa,cex=0.7,main="2 Break points",xlab="Age",ylab="lpsa") 
 lines(sort(data$age),pred3,type="p",col="blue",pch=15,cex=0.7)
plot(data$age,data$lpsa,cex=0.7,main="100 Break points",xlab="Age",ylab="lpsa")
 lines(sort(data$age),pred4,type="p",col="purple",pch=15,cex=0.7)

Adesso proviamo ad usare lcavol, nochè la variabile indipendente con la correlazione più alta rispetto alla nostra variabile risposta:

cut<-cut(sort(data$lcavol),4)
mcut<-lm(lpsa~cut,data)
pred<-predict(mcut,data)

cut2<-cut(sort(data$lcavol),10)
mcut2<-lm(lpsa~cut2,data)
pred2<-predict(mcut2,data)

cut3<-cut(sort(data$lcavol),2)
mcut3<-lm(lpsa~cut3,data)
pred3<-predict(mcut3,data)

cut4<-cut(sort(data$lcavol),100)
mcut4<-lm(lpsa~cut4,data)
pred4<-predict(mcut4,data)

par(mfrow=c(2,2))
plot(data$lcavol,data$lpsa,cex=0.7,main="4 Break points",xlab="lcavol",ylab="lpsa")
 lines(sort(data$lcavol),pred,type="p",col="red",pch=15,cex=0.7)
plot(data$lcavol,data$lpsa,cex=0.7,main="10 Break points",xlab="lcavol",ylab="lpsa") 
 lines(sort(data$lcavol),pred2,type="p",col="green",pch=15,cex=0.7)
plot(data$lcavol,data$lpsa,cex=0.7,main="2 Break points",xlab="lcavol",ylab="lpsa") 
 lines(sort(data$lcavol),pred3,type="p",col="blue",pch=15,cex=0.7)
plot(data$lcavol,data$lpsa,cex=0.7,main="100 Break points",xlab="lcavol",ylab="lpsa")
 lines(sort(data$lcavol),pred4,type="p",col="purple",pch=15,cex=0.7)

Le step function, pur potendo cambare il numero di break points, rimangono molto grezze e poco utili per la stima di valori, ma sono sicuramente un ottimo metodo per un’analisi primitiva dei dati: infatti viene totalmente colta da tutti i modelli la tendenza crescente dei dati e nel caso dei modelli con più break points (tranne il modello con 100 punti fatto più per curiosità che per utilità) viene colta anche la differenza di velocità di crescita tra un intervallo e l’altro.

Forest plot di una regressione di tipo “Lasso”

Proviamo a creare un forest plot per un modello di tipo Lasso, come fatto in precedenza con i modello lineare e polinomiale.

library(glmnet)
## Loaded glmnet 4.1-6
library(ggplot2)
set.seed(1)
x<-data.matrix(data[,-9])
nomivar <- c("lcavol","lweight","age","lbph","svi","gleason","pgg45")
cv.lambda<-cv.glmnet(x,data$lpsa,alpha=1,nfolds=30)
lambda<- cv.lambda$lambda.min

l<-glmnet(x,data$lpsa,alpha=1,lambda=lambda)
coeff<-l[["beta"]]@x

df<-data.frame(var=nomivar,coef=coeff)
ggplot(data=df, aes(y=var, x=coef)) +
  geom_point()  +
  labs( x='Coefficients', y = 'Predictor') +
  geom_vline(xintercept=0, color='black', linetype='dashed', alpha=.5) +
  theme_minimal()+
  expand_limits(x = c(-0.5,1))

A differenza dei modelli precedenti, il forest plot appena creato mostra chiaramente il funzionamento di una regressione Lasso: le covariate non significative secondo il modello vengono “shrinkate” a 0 in maniera decisa (guardiamo age, gleason ,pgg45 e lbph) rispetto alle altre che possono avere un ruolo più importante all’interno del modello (lweight, lcavol e svi).

N.B.

Ho provato ad inserire all’interno di questo forest plot i vari standard error dei coefficient, ma nel cercare una soluzione “manuale” a questo problema sono incappato in libri e paper nella quale si discuteva dell’effettiva utilità degli standard error in modelli di regressione penalizzanti come il Lasso. Rimane però che al momento non abbia ancora trovato un modo per implementare questa informazione all’interno dell’analisi, ma cercherò di trovare una soluzione per i prossimi laboratori.

Regressione locale

Infine, cerchiamo di stimare delle regressini locali per i nostri dati con un modello semplificato, composto solo da una variabile risposta, lpsa, e una covariata, lcavol:

rl5<-loess(data$lpsa~data$lcavol,data,span=0.5)
pred5<-predict(rl5,sort(data$lcavol))
plot(data$lcavol,data$lpsa,ylab="lpsa",xlab="lcavol",main="Regressioni locali con diversi parametri",cex=0.8)
legend("bottomright",legend=c("span=0.5","span=0.1","span=0.9"),col=c("blue","red","green"),
       lty=1,cex=0.8,box.lty=0,lwd=2)
lines(sort(data$lcavol),pred5,lwd=2,col="blue")

rl1<-loess(data$lpsa~data$lcavol,data,span=0.1)
pred1<-predict(rl1,sort(data$lcavol))
lines(sort(data$lcavol),pred1,col="red",lwd=2)

rl9<-loess(data$lpsa~data$lcavol,data,span=0.9)
pred9<-predict(rl9,sort(data$lcavol))
lines(sort(data$lcavol),pred9,col="green",lwd=2)

Usando diversi valori per il parametro span possiamo vedere l’effetto che questo provoca alla nostra regressione: più il valore è vicino a 0 e più la regressione tenderà a interpolare più dati possibili creando una curva spigolosa che varia in maniera selvaggia in base ai punti; per valori vicini a 1 ( ma già a 0.5 è apprezzabile) la regressione risulta più “dolce” e meno influenzata dalla posizione dei singoli punti, ma piuttosto da una “media” di questi ultimi.