Esempio per usare Rmarkdown

Innanzitutto dovrete avere installato il pacchetto Rmarkdown.

Per questo esempio rifaremo quanto svolto nel laboratorio in classe con i dati tratti da mcycle. Guardate con attenzione l’output di questo programma e intuirete l’uso di alcune convenzioni per mettere i titoli ai paragrafi, per metter il grassetto o l’italico e soprattutto per gestire il codice R inserendo i chunk (ovvero i pezzi di programma)

Si noti che dopo avere introdotto il testo e i chunk, sarà necessario in Rstudio cliccare sul tasto knit che compare nella finestra in cui si scrive il codice Rmarkdown (che verrà poi salvato automaticamente)

Prendiamo i dati

La prima cosa (si noti che saltando una linea si andra a capo) sarà di assicurarsi di avere caricato i pacchetti (packages) che sono necessari usando la funzione install.packages("nome package") (si possono anche installare preventivamente mediante il menù in Rstudio.

I dati sono contenuti nella libreria MASS (si noti come faccio a ottenere il carattere che si usa nello script R). Per l’esercizio occorre quindi i dati con i comandi che vedete

# quanto sopra apre un chunk di comandi R. Il nome di questo chunk è "dati" 
# (può essere anche omesso). 
# Fra le parentesi di apertura del chunk si possono aggiungere parametri, ad esempio
# mettesse echo=off allora nel documento non verranno riportati i comandi R. 

library(MASS)
data(mcycle)

# possiamo anche far un primo scatterplot dei dati. Che comparirà nel report

plot(mcycle$times, mcycle$accel, ylab="accellerazione", xlab="tempi")

Generiamo il training e il test set

Il primo passo sarà generare:

il training set
il test set

dividendo casualmente il nostro insieme di dati. (si noti come si è gestito l’elenco puntato)

n<-length(mcycle$accel)
set.seed(4231)
train<-sample(n,round(0.75*n))
train_set<-mcycle[train,]
test_set<-mcycle[-train,]

Adattiamo un modello lineare

Iniziamo adattando un modello lineare, dal plot visto non ci aspettermo che vada granchè bene. Il modello verrà stimato solo sui dati del training set.

La funzione che utilizzeremo per il modello lineare sarà lm()

mod1<-lm(accel~times, data=train_set)

Ora dovremo calcolare l’errore di adattamento e quello di previsione (o generalizzazione). Per fare questo calcoleremo le previsioni prima sul training set e poi sul test set e quindi calcoleremo le due devianze medie con la formula (si guardi questo primo esempio per gestire le formule)

\(1/n \sum_{i=1}^n (y_i-\hat{y}_i)^2\)

fit1ad<-predict(mod1, newdata = train_set) #previsioni sul training set
fit1pr<-predict(mod1, newdata = test_set) ##previsioni sul test set

# ora calcoliamo le devianze medie nei due casi
errad1<-mean((fit1ad-train_set$accel)^2)
errpr1<-mean((fit1pr-test_set$accel)^2)

Modelli più complessi

Proviamo ora modelli più complessi. Prima un modello cubico e calcoliamo ancora errore di adattamento e di previsione

mod3<-lm(accel~poly(times, 3), data=train_set)
fit3ad<-predict(mod3, newdata = train_set)
fit3pr<-predict(mod3, newdata = test_set)
errad3<-mean((fit3ad-train_set$accel)^2)
errpr3<-mean((fit3pr-test_set$accel)^2)

Proviamo un modello di grado 6 e calcoliamo errore di adattamento e di previsione

mod6<-lm(accel~poly(times, 6), data=train_set)
fit6ad<-predict(mod6, newdata = train_set)
fit6pr<-predict(mod6, newdata = test_set)
errad6<-mean((fit6ad-train_set$accel)^2)
errpr6<-mean((fit6pr-test_set$accel)^2)

Proviamo un modello di grado 15 e poi di grado 20 e calcoliamo errore di adattamento e di previsione

#prima grado 15
mod15<-lm(accel~poly(times, 15), data=train_set)
fit15ad<-predict(mod15, newdata = train_set)
fit15pr<-predict(mod15, newdata = test_set)
errad15<-mean((fit15ad-train_set$accel)^2)
errpr15<-mean((fit15pr-test_set$accel)^2)

# e poi un modello di grado 20 e calcoliamo errore di adattamento e di previsione

mod20<-lm(accel~poly(times, 20), data=train_set)
fit20ad<-predict(mod20, newdata = train_set)
fit20pr<-predict(mod20, newdata = test_set)
errad20<-mean((fit20ad-train_set$accel)^2)
errpr20<-mean((fit20pr-test_set$accel)^2)

Per capire cosa accade al crescere della complessità otteniamo un grafico con i diversi modelli polinomiali

plot(accel~times, cex=1, data=train_set)
points(accel~times, data=test_set, cex=2, pch=20, col=2)
abline(mod1, lwd=2, col=3)
legend(45,-75, c("test set", "training set"), pch=c(1,20), col=c(1,2))
lines(seq(0,60,.05),predict(mod3,data.frame(times=seq(0,60,.05))), lwd=2, col=3)
lines(seq(0,60,.05),predict(mod6,data.frame(times=seq(0,60,.05))), lwd=2, col=4)
lines(seq(0,60,.05),predict(mod15,data.frame(times=seq(0,60,.05))), lwd=2, col=5)
lines(seq(0,60,.05),predict(mod20,data.frame(times=seq(0,60,.05))), lwd=2, col=6)

Ora possiamo vedere quale andamento hanno gli errori al crescere del grado

ad<-c(errad1,errad3,errad6,errad15,errad20)
pr<-c(errpr1,errpr3,errpr6,errpr15, errpr20)
plot(c(1,3,6,15,20), pr, xlab="grado del polinomio", ylab="errori", pch=20, ylim=c(0,2800))
points(c(1,3,6,15,20), ad,  pch=20, col=2)
legend(11,2000, c("errore di adattamento", "errore di previsione"),pch=c(20,20), col=c(2,1))

Utilizzo di un loop con ciclo for

Ora proviamo a impostare tutto con un ciclo for. Salveremo i valori in due vettori e faremo variare il grado da 1 a 20

k=22
errad<-rep(1,k)
errpr<-rep(1,k)

for(i in 1:k){
  mod<-lm(accel~poly(times, i), data=train_set)
  fitad<-predict(mod, newdata = train_set)
  fitpr<-predict(mod, newdata = test_set)
  errad[i]<-mean((fitad-train_set$accel)^2)
  errpr[i]<-mean((fitpr-test_set$accel)^2)
}

Infine proviamo a ottenere il grafico dei due eroori al crescere della complessità del modello

plot(1:k, errad[1:k], pch=20, col=2, xlab="grado del polinomio", ylab="errori", ylim=c(200,2500))
points(1:k, errpr[1:k], pch=20)
legend(13,1500, c("errore di adattamento", "errore di previsione"),pch=c(20,20), col=c(2,1))

Esempio