Punto 1: generazione dei dati da una spline cubica

Avendo a disposizione il vettore dei coefficienti, i nodi e un vettore di x (generato da un’uniforme(-1,1) geenriamo delle y da una spline cubica a basi troncate:

library(splines)
library(data.table)
library(ggplot2)
library(dplyr)
set.seed(77)
x<-runif(150,min=-1,max=1)
e<-rnorm(150,0,1)

b<-c(1,0.7,-0.1,5,-1.5,0.1,-0.3)
knots<-c(-0.7,0,0.7)

y<-c()

for(i in 1:150){
  y[i]<-b[1]+b[2]*x[i]+b[3]*(x[i]^2)+b[4]*(x[i]^3)
  
  if ((x[i]-knots[1])>0){
    y[i]<-y[i]+b[5]*(x[i]-knots[1])^3
  }
  if ((x[i]-knots[2])>0){
    y[i]<-y[i]+b[6]*(x[i]-knots[2])^3
  }
  if ((x[i]-knots[3])>0){
    y[i]<-y[i]+b[7]*(x[i]-knots[3])^3
  }
  else{
    y[i]<-y[i]
  }
 y[i]<-y[i]+e[i]
}

plot(x,y,main="Dati generati da una spline cubica con 3 nodi fissati")

Al vettore delle y aggiungiamo un valore e per simulare un effetto dato da un errore stocastico, altrimenti i punti si disporebbero perfetamente come la spline dalla quale sono generati.

plot(x,(y-e),main="Dati senza errore stocastico",ylab="y(senza errore stocastico)")

Punto 2: stima di una spline sui dati generati

Ora proviamo a strimare una prima spline cubica dove specifichiamo i nodi e una seconda dove verranno specificati i gradi di libertà. Nel caso della seconda spline specifichiamo 6 gradi di libertà e non 7 (K+4) perchè il parametro df richiede il numero di gdl-1 se è presente l’intercetta e nel nostro caso è così.

d<-data.frame(x,y)
d<-arrange(d,x)
s1<-lm(d$y~bs(d$x,knots=c(-0.7,0,0.7)),d)
plot(x,y,cex=0.7,main="Spline cubiche a confronto")
pred<-predict(s1,d)
lines(d$x,pred,col="red",lwd=2)
s2<-lm(d$y~bs(d$x,df=6),d)
pred2<-predict(s2,d)
lines(d$x,pred2,col="blue",lwd=2)
legend("bottomright",legend=c("Nodi fissati","Gdl fissati"),col=c("red","blue"),lwd=2,box.lty=0)

Non si notano grandi differenze tra le due curve, se non una leggerissima tendenza a “divergere” nelle code della spline con i nodi fissati rispetto alla spline con i gdl fissati.

Punto 3: stima di spline naturali e di lisciamento

Stimiamo ora spline naturali e di lisciamento provando diversi valori per il parametro di tuning \(\lambda\)

plot(d,cex=0.7,main="Spline naturali cubiche con diversi Gdl")
natural1<-lm(d$y~ns(d$x,df=1),d)
lines(d$x,predict(natural1,d),col="green",lwd=2)

natural2<-lm(d$y~ns(d$x,df=7),d)
lines(d$x,predict(natural2,d),col="blue",lwd=2)

natural3<-lm(d$y~ns(d$x,df=3),d)
lines(d$x,predict(natural3,d),col="orange",lwd=2)

legend("bottomright",legend=c("1 Gdl","3 Gdl","7 Gdl"),col=c("green","orange","blue"),lwd=2,box.lty=0)

Più gradi di libertà si fissano e più la curva è “flessibile”, con un solo Gdl la curva è una linea mentre già con 3 Gdl assomiglia ad una polinomiale in grado di interpolare i dati.

plot(d,cex=0.7,main=expression(paste("Spline di lisciamento con diversi valori ",lambda)))
smooth1<-smooth.spline(d$x,d$y,lambda=1)
lines(smooth1,col="red",lwd=2)

smooth2<-smooth.spline(d$x,d$y,lambda=0.5)
lines(smooth2,col="green",lwd=2)

smooth3<-smooth.spline(d$x,d$y,lambda=0.005)
lines(smooth3,col="purple",lwd=2)

smooth4<-smooth.spline(d$x,d$y,lambda=1000)
lines(smooth4,col="blue",lwd=2)

legend("bottomright",legend=c(expression(paste(lambda,"=1")),expression(paste(lambda,"=0.5")),expression(paste(lambda,"=0.005")),expression(paste(lambda,"=1000"))),col=c("red","green","purple","blue"),lwd=2,box.lty=0)

Si può vedere come al variare del parametro \(\lambda\) cambiano le forme delle curve: più \(\lambda\) è grande e più la spline di lisciamento viene regolarizzata (esempio estremo con \(\lambda\)=1000 dove la curva diventa una linea retta)

Punto 4: da basi troncate a B-splines

Ristimiamo tutte le varie spline precedenti, ma stavolta non tramite basi troncate ma tramite B-splines:

basi<-bs(x, knots = knots, intercept = TRUE)

y2<- basi %*% b
plot(x,y2,main="Dati generati da B-splines senza errore stocastico",ylab="y")

y2<-y2+e
plot(x,y2,main="Dati generati da B-splines con errore stocastico",ylab="y")

Anche in questo caso aggiungiamo il vettore e per simulare l’errore stocastico all’interno dei nostri dati

d2<-data.frame(x,y2)
d2<-arrange(d2,x)
plot(d2,cex=0.7,main="Spline cubiche a confronto",ylab="y")
s1.2<-lm(d2$y2~bs(d2$x,knots=c(-0.7,0,0.7)),d2)
pred<-predict(s1.2,d2)
lines(d2$x,pred,col="red",lwd=2)

s2.2<-lm(d2$y~bs(d2$x,df=6),d2)
pred2<-predict(s2.2,d2)
lines(d2$x,pred2,col="blue",lwd=2)

legend("topright",legend=c("Nodi fissati","Gdl fissati"),col=c("red","blue"),lwd=2,box.lty=0)

plot(d2,cex=0.7,main="Spline naturali cubiche con diversi Gdl",ylab="y")
natural1<-lm(d2$y~ns(d2$x,df=1),d2)
lines(d2$x,predict(natural1,d2),col="green",lwd=2)

natural2<-lm(d2$y~ns(d2$x,df=5),d2)
lines(d2$x,predict(natural2,d2),col="blue",lwd=2)

natural3<-lm(d2$y~ns(d2$x,df=3),d2)
lines(d2$x,predict(natural3,d2),col="orange",lwd=2)

legend("topright",legend=c("1 Gdl","3 Gdl","7 Gdl"),col=c("green","orange","blue"),lwd=2,box.lty=0)

plot(d2,cex=0.7,main=expression(paste("Spline di lisciamento con diversi valori ",  lambda)),ylab="y")
smooth1<-smooth.spline(d2$x,d2$y,lambda=1)
lines(smooth1,col="red",lwd=2)

smooth2<-smooth.spline(d2$x,d2$y,lambda=0.5)
lines(smooth2,col="green",lwd=2)

smooth3<-smooth.spline(d2$x,d2$y,lambda=0.005)
lines(smooth3,col="purple",lwd=2)

smooth4<-smooth.spline(d2$x,d2$y,lambda=1000)
lines(smooth4,col="blue",lwd=2)

legend("topright",legend=c(expression(paste(lambda,"=1")),expression(paste(lambda,"=0.5")),expression(paste(lambda,"=0.005")),expression(paste(lambda,"=1000"))),col=c("red","green","purple","blue"),lwd=2,box.lty=0)

Punto 5: stime sui dati “Spline_data”

data<-read.csv("splines_data.csv",header=TRUE,row.names=1)
par(mfrow=c(2,2))
plot(data)
hist(data$X)
hist(data$Y)

  
m<-lm(data$Y~ data$X,data)
plot(data,cex=0.7)
abline(m,col="red",lwd=2)

Ovviamente i dati non sono per nulla interplabili da una regressione lineare, usiamola come confronto per apprezzare meglio le stime che si faranno successivamente con le spline

#splines cubiche
plot(data,cex=0.7,main="Spline cubiche con diversi gradi di libertà")

bs1<-lm(data$Y~bs(data$X,df=3),data)
bs2<-lm(data$Y~bs(data$X,df=5),data)
bs10<-lm(data$Y~bs(data$X,df=10),data)
bs20<-lm(data$Y~bs(data$X,df=20),data)

lines(data$X,predict(bs1,data),col="purple",lwd=2)
lines(data$X,predict(bs2,data),col="red",lwd=2)
lines(data$X,predict(bs10,data),col="blue",lwd=2)
lines(data$X,predict(bs20,data),col="green",lwd=2)

legend("bottomright",legend=c("3 Gdl","5 Gdl","10 Gdl","20 Gdl"),col=c("purple","red","blue","green"),lwd=2,box.lty=0)

#natural splines
plot(data,cex=0.7,main="Natural splines con diversi gradi di libertà")
ns1<-lm(data$Y~ns(data$X,df=3),data)
ns2<-lm(data$Y~ns(data$X,df=5),data)
ns10<-lm(data$Y~ns(data$X,df=10),data)
ns20<-lm(data$Y~ns(data$X,df=20),data)

lines(data$X,predict(ns1,data),col="purple",lwd=2)
lines(data$X,predict(ns2,data),col="red",lwd=2)
lines(data$X,predict(ns10,data),col="blue",lwd=2)
lines(data$X,predict(ns20,data),col="green",lwd=2)

legend("bottomright",legend=c("3 Gdl","5 Gdl","10 Gdl","20 Gdl"),col=c("purple","red","blue","green"),lwd=2,box.lty=0)

#smooth spline
plot(data,cex=0.7,main=expression(paste("Spline di lisciamento con diversi valori ",  lambda)))

sm1<-smooth.spline(data$X,data$Y,lambda=1)
sm0.005<-smooth.spline(data$X,data$Y,lambda=0.005)
smlim<-smooth.spline(data$X,data$Y,lambda=0.00005)
sm1000<-smooth.spline(data$X,data$Y,lambda=1000)


lines(sm1,col="red",lwd=2)
lines(sm0.005,col="blue",lwd=2)
lines(smlim,col="darkgreen",lwd=2)
lines(sm1000,col="purple",lwd=2)


legend("bottomright",legend=c(expression(paste(lambda,"=1")),expression(paste(lambda,"=0.005")),expression(paste(lambda,"=0.0005")),expression(paste(lambda,"=1000"))),col=c("red","blue","darkgreen","purple"),lwd=2,box.lty=0)

Punto 6: GAM sui dati trees

library(mgcv)
## Caricamento del pacchetto richiesto: nlme
## 
## Caricamento pacchetto: 'nlme'
## Il seguente oggetto è mascherato da 'package:dplyr':
## 
##     collapse
## This is mgcv 1.8-41. For overview type 'help("mgcv-package")'.
data("trees")

summary(trees)
##      Girth           Height       Volume     
##  Min.   : 8.30   Min.   :63   Min.   :10.20  
##  1st Qu.:11.05   1st Qu.:72   1st Qu.:19.40  
##  Median :12.90   Median :76   Median :24.20  
##  Mean   :13.25   Mean   :76   Mean   :30.17  
##  3rd Qu.:15.25   3rd Qu.:80   3rd Qu.:37.30  
##  Max.   :20.60   Max.   :87   Max.   :77.00
plot(trees)

hist(trees$Volume)

hist(log(trees$Volume))

trees$Volume<-log(trees$Volume)

Dopo aver ddato un’occhiata ai dati ed aver traformato la variabileVolume, passiamo alle stime dei modelli GAM:

#parametro sp indica il lambda delle funzioni di lisciamento
gam1<-gam(trees$Volume~s(trees$Girth,sp=0.45),data=trees)
gam2<-gam(trees$Volume~s(sort(trees$Height),sp=0.45),data=trees)
gam3<-gam(trees$Volume~s(trees$Girth,sp=0.45)+ s(sort(trees$Height),sp=0.45),data=trees)

pred1<-predict(gam1)
pred2<-predict(gam2)
pred3<-predict(gam3)

plot(trees$Girth,trees$Volume,cex=0.7,xlab="Girth",ylab="log(Volume)")
lines(trees$Girth,pred1,col="red",lwd=2)

plot(trees$Height,trees$Volume,cex=0.7,ylab="log(Volume)")
lines(sort(trees$Height),pred2,col="blue",lwd=2)

plot(trees$Girth,trees$Volume,cex=0.7,xlab="Girth",ylab="log(Volume)")
lines(trees$Girth,pred3,col="green",lwd=2)

Confromtiamo tra loro i modelli attraverso un test Anova:

anova(gam1,gam3)
anova(gam2,gam3)

Il test ci suggerisce, per entrambi i casi, che tra il modello semplice con una covariata e quello “più” complesso non ci sia molta differenza, o almeno, tale da giustificare l’uso di un modello più complicato a beneficio di una magggior efficacia

summary(gam3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## trees$Volume ~ s(trees$Girth, sp = 0.45) + s(sort(trees$Height), 
##     sp = 0.45)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.27273    0.02071     158   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                         edf Ref.df      F  p-value    
## s(trees$Girth)        1.822  2.257 17.406 1.12e-05 ***
## s(sort(trees$Height)) 2.185  2.701  1.131    0.287    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.952   Deviance explained = 95.8%
## GCV = 0.015856  Scale est. = 0.013295  n = 31

Dal summary del modello gam3 possiamo notare come gli effetti parametrici siano predominanti rispetto a quelli non parametrici, questro ci indica l’assenza di relazioni non lineari all’interno del modello, a riprova di ciò possiamo confrontarlo con un modello lineare per vederne l’estrema somiglianza

m<-lm(Volume~ Girth+Height,data=trees)
predm<-predict(m)
par(mfrow=c(1,2))
plot(trees$Girth,trees$Volume,cex=0.7,xlab="Girth",ylab="log(Volume)",main="Modello lineare")
lines(trees$Girth,predm)
plot(trees$Girth,trees$Volume,cex=0.7,xlab="Girth",ylab="log(Volume)",main="Modello GAM")
lines(trees$Girth,pred3,col="green",lwd=2)