---
title: "Esempio"
output: html_document
date: "2024-10-07"
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
# 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
```{r dati}
# 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)
```{r traintest }
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()```
```{r linmod}
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$
```{r errlinmod}
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
```{r cubic}
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
```{r poli6}
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
```{r gradoalto}
#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
```{r grafici}
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
```{r errori}
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
```{r ciclo}
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
```{r graficoerr}
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))
```