R dispone di una enorme numero di funzionalità aggiuntive al nucleo fondamentale del linguaggio, che sono organizzate sotto forma di package. Ecco un elenco di funzioni per manipolare i package.
library(): visualizza l'elenco dei package disponibili nel sistema;library(package): carica in memoria il package specificato;help(...., package): visualizza un breve messaggio di informazioni sul package indicato;help(topic, ... ,package): visualizza messaggi di informazione sulla funzione topic all'interno del package specificato. Se il package è stato caricato in memoria, l'opzione package è superflua.Nel package rpart è implementato un algoritmo per l'induzione di alberi di classificazione e di regressione.
rpart(formula, data, ...): determina un albero di classificazione / regressione, operando sul data frame indicato nel parametro data, secondo la formula specificata. La formula è un costrutto del tipo vardip ~ var1 + .... + varn dove vardip, var1, ... varn sono elementi del data frame di input, e indica che vardip è la variabile classe che si vuole predirre in base alle variabili esplicative var1, ..., varn. Il risultato è una lista di classe rpart. È possibile modificare il funzionamento di rpart con alcuni parametri:
minsplit: normalmente uguale a 20. Se un nodo ha meno istanze associate di minsplit, non verrà mai diviso;minbucket: normalmente uguale a minsplit/3, è il minimo numero di istanze che può essere associato ad una foglia.predict(object,newdata=list()): applica il modello in object ai dati in newdata.plot(x): se x è un oggetto di classe rpart, visualizza l'albero di classificazione corrispondente, ma senza etichette.text(x): se x è un oggetto di classe rpart, visualizza le etichette dell'albero di classificazione.Ad esempio, consideriamo la seguente sessione di lavoro sul dataset iris:
> model=rpart(Species ~ .,iris)
> model
n= 150
node), split, n, loss, yval, (yprob)
* denotes terminal node
1) root 150 100 setosa (0.33333333 0.33333333 0.33333333)
2) Petal.Length< 2.45 50 0 setosa (1.00000000 0.00000000 0.00000000) *
3) Petal.Length>=2.45 100 50 versicolor (0.00000000 0.50000000 0.50000000)
6) Petal.Width< 1.75 54 5 versicolor (0.00000000 0.90740741 0.09259259) *
7) Petal.Width>=1.75 46 1 virginica (0.00000000 0.02173913 0.97826087) *
> plot(model)
> text(model)
> model=rpart(Species ~ .,iris,minbucket=2)
> model
n= 150
node), split, n, loss, yval, (yprob)
* denotes terminal node
1) root 150 100 setosa (0.33333333 0.33333333 0.33333333)
2) Petal.Length< 2.45 50 0 setosa (1.00000000 0.00000000 0.00000000) *
3) Petal.Length>=2.45 100 50 versicolor (0.00000000 0.50000000 0.50000000)
6) Petal.Width< 1.75 54 5 versicolor (0.00000000 0.90740741 0.09259259)
12) Petal.Length< 4.95 48 1 versicolor (0.00000000 0.97916667 0.02083333) *
13) Petal.Length>=4.95 6 2 virginica (0.00000000 0.33333333 0.66666667) *
7) Petal.Width>=1.75 46 1 virginica (0.00000000 0.02173913 0.97826087) *
Nell'output della variabile model abbiamo una serie di righe, ognuna della forma "node), split, n, loss, yval, (yprob)", che danno informazioni sull'albero generato. L'indentazione evidenzia la relazione genitore-figlio tra i nodi, mentre i dati della singola riga hanno questo significato:
node: numero del nodo. Ogni nodo x, tranne il nodo 1, ha come padre il nodo floor(x/2), dove floor(x) è la parte intera inferiore di x (ad esempio, il nodo 13 ha come padre il nodo floor(6.5)=6);split: la condizione per cui si arriva a quel nodo a partire dal genitore;n: il numero di istanze che fanno capo al nodo;yprob: la distribuzione di probabilità delle classi per le istanze associate al nodo;yval: la classe più numerosa nelle istanze associate al nodo;loss: il numero di istanze la cui classe non coincide con yval.Le foglie sono inoltre etichettate con un asterisco. Per sapere il numero di errori compiuti nell'insieme di addestramento, basta sommare il valore di loss per tutti i nodi foglia.
> iris2=data.frame(Sepal.Length=c(6.1,5.1),Sepal.Width=c(3.3,2.3),Petal.Length=c(5.4,3.9),Petal.Width=c(2,1.1));
> plot(iris$Petal.Length,iris$Petal.Width,col=as.integer(iris$Species),pch=as.integer(iris$Species))
> points(iris2$Petal.Length,iris2$Petal.Width,col="blue",pch="X")
> predict(model,iris2)
setosa versicolor virginica
[1,] 0 0.02173913 0.97826087
[2,] 0 0.97916667 0.02083333
Dal diagramma si può vedere come la prima istanza sia nel mezzo del territorio dell'iris virginica, mentre la seconda sia nel territorio dell'iris versicolor. Applicando il metodo predict si vede che queste considerazioni visuali trovano conferma dal risultato del classificatore: ad ogni riga corrisponde una distribuzione di probabilità per la corrispondente istanza, ottenuta semplicemente copiando la distribuzione di probabilità della foglia del modello che corrisponde all'istanza.
Il package class implementa l'algoritmo k-nearest neighbour. La distanza utilizzata è quella euclidea, e funziona solo con attributi numerici. Le funzioni da usare sono le seguenti:
knn(train, test, cl, k=1,l=0, prob=FALSE, use.all=TRUE): applica l'algoritmo k-nn hai dati di test sulla base dei dati di addestramento. Vediamo in dettaglio i vari parametri:
train è un data-frame con l'insieme di addestramento (senza l'attributo classe);test è un data-frame con l'insieme di test (senza attributo classe) o un singolo vettore nel caso si voglia applicare l'algoritmo su una sola istanza;k è l'ovvio parametro k dell'algoritmo;cl è un vettore con l'attributo classe dei dati di addestramento. Accetta solo attributi numerici;l contiene il numero minimo di voti per considerare significativo un risultato, altrimenti viene prodotto il risultato speciale doubt;prob se pari a TRUE aggiunge nell'attributo prob del risultato le distribuzioni di probabilità delle varie classi;use.all determina cosa fare nel caso di più istanze alla stessa distanza; se è pari a TRUE tutti vengono considerato nel determinare il valore della nuova istanza, anche a costo di superare il valore k; se è pari a FALSE, se necessario viene estratto casualmente un campione di essi che porti esattamente a k il numero di vicini considerati.?knn.
knn.cv(train, cl, k=1,l=0, prob=FALSE, use.all=TRUE): use il metodo del LOO croos-validation per testare la bontà del metodo k-nearest neighbour. Il risultato è il valore di classe ottenuto per ogni istanza dell'insieme di addestramento sulla base di tutte le altre.knn.cv(iris[-5],iris[[5]])
Confrontare la precisione di knn e rpart sul data-set dell'iris.
Esistono varie funzioni in R che consentono di effettuare la regressione lineare, semplice o multipla. Il comando di base è il seguente:
lm(formula, data=df): esegue il calcolo della regressione lineare, sul set di dati indicato nel parametro data, secondo la formula specificata. La formula è un costrutto del tipo vardip ~ var1 + .... + varn dove vardip, var1, ... varn sono elementi del data frame di input, e indica la dipendenza lineare di vardip dalle variabili esplicative var1, ..., varn.
Il risultato è l'elenco dei parametri della retta di regressione, ma può essere salvato per essere ulteriormente analizzato con altre funzioni.
> florida <- read.table("florida.txt");
> lm(BUCHANAN ~ BUSH, data=florida)
Call:
lm(formula = BUCHANAN ~ BUSH, data = florida)
Coefficients:
(Intercept) BUSH
45.289861 0.004917
In pratica, il risultato di lm è una lista con varie componenti che contengono, oltre ai parametri della regressione, altre informazioni aggiuntive. Normalmente, però, R ci visualizza soltanto i parametri di regressione. Se vogliamo visualizzare la lista completa, possiamo usare il comando unclass:
> florida <- read.table("florida.txt");
> unclass(lm(BUCHANAN ~ BUSH, data=florida))
$coefficients
(Intercept) BUSH
45.289861271 0.004916828
$residuals
1 2 3 4 5 6
49.2331397 0.1267330 12.7386511 -6.9046518 -41.6347069 -127.9402315
7 8 9 10 11 12
.... omissis ....
Il risultato di lm viene tipicamente salvato e usato come input ad una delle seguenti funzioni, che consentono di esplorare ulteriormente il modello trovato dalla regressione lineare:
summary(lm): visualizza informazioni aggiuntive sul risultato di una regressione lineare. Il significato dei vari campi di output si può trovare in libri di statistica.
> model <- lm(BUCHANAN ~ BUSH, data=florida)
> summary(model)
Call:
lm(formula = BUCHANAN ~ BUSH, data = florida)
Residuals:
Min 1Q Median 3Q Max
-907.50 -46.10 -29.19 12.26 2610.19
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.529e+01 5.448e+01 0.831 0.409
BUSH 4.917e-03 7.644e-04 6.432 1.73e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 353.9 on 65 degrees of freedom
Multiple R-Squared: 0.3889, Adjusted R-squared: 0.3795
F-statistic: 41.37 on 1 and 65 DF, p-value: 1.727e-08
coef(lm): restiuisce un vettore con i coefficienti di regressione del modello lm.
> coef(model) (Intercept) BUSH 45.289861271 0.004916828
residuals(lm): restituisce un vettore con i residui (differenza tra il valore velore dell'attributo classe e il valore predetto) di regressione del modello lm.
> residuals(model)
1 2 3 4 5 6
49.2331397 0.1267330 12.7386511 -6.9046518 -41.6347069 -127.9402315
7 8 9 10 11 12
30.5840916 -37.4389960 78.4640035 -64.5428507 -220.3941165 -10.1979647
13 14 15 16 17 18
-907.4952582 -30.2158817 -29.5554635 -143.0509136 99.6390988 -24.2812301
....omissis...
fitted(lm): restituisce un vettore con i valori calcolati dal modello lm, sullo stesso insieme di dati di addestramento.
> residuals(model)
1 2 3 4 5 6
49.2331397 0.1267330 12.7386511 -6.9046518 -41.6347069 -127.9402315
7 8 9 10 11 12
30.5840916 -37.4389960 78.4640035 -64.5428507 -220.3941165 -10.1979647
13 14 15 16 17 18
-907.4952582 -30.2158817 -29.5554635 -143.0509136 99.6390988 -24.2812301
...omissis...
predict(lm,newdf): utilizza il modello lineare in lm per predirre i nuovi dati presenti nel dataframe newdf. Quest'ultimo deve avere tutti i campi utilizzati come variabili indipendenti in lm.
> predict(model,data.frame(BUSH=c(150000,120000)))
1 2
782.8141 635.3092
Alcuni set di dati su cui poter fare esperimenti sono:
Usando la funzione di regressione lineare, esaminiamo i risultati delle elezioni presidenziali USA del 2000. Vi ricordo che queste elezioni hanno suscitato molto scalpore perché in una delle circoscrizioni della Florida (Miami-Dade), il posto che normalmente era occupato da Gore nella scheda elettorale era invece occupato da Buchanan. A detta dei sostenitori di Gore, questo potrebbe aver causato un certo numero di errori nella votazione, che hanno consentito a Bush di vincere le elezioni in Florida, e quindi in tutti gli USA. Vediamo se questa ipotesi è plausibile.
Esaminiamo graficamente la relazione tra i voti di Bush e quelli di Buchanan nelle varie contee.
> florida = read.table("florida.txt")
> attach(florida)
> plot(BUSH, BUCHANAN)
Vediamo che ci sono due contee (una in particolare) che si comportano da outliers. Identifichiamo queste contee con il comando identify.
> identify(BUSH, BUCHANAN, n=2)
[1] 13 50
> florida[c(13,50),]
County GORE BUSH BUCHANAN NADER BROWN HAGELIN HARRIS MCREYNOLDS MOOREHEAD PHILLIPS Total
13 DADE 328702 289456 561 5355 759 119 88 36 124 69 625269
50 PALM BEACH 268945 152846 3407 5564 743 143 45 302 103 188 432286
Una delle due contee "strane" è proprio quella di Dade. Supponiamo allora che effettivamente ci sia stato un errore consistente nella procedura di votazione, e proviamo a calcolare il numero di voti corretto per Buchanan con la regressione lineare.
> model = lm(BUCHANAN ~ BUSH, data=florida[-50,])
> abline(coef(model)[1],coef(model)[2],col="red")
> predict(model, data.frame(BUSH=BUSH[50]))
1
597.7677
> BUCHANAN[50]
[1] 3407
Dunque, secondo i nostri conti, Buchanan avrebbe dovuto prendere circa 598 voti nella contea di Dade, e non 3407. La differenza è di 2809. Se questi fossero veramente voti "perduti" da Gore, siccome la differenza tra i voti di Bush e di Gore in Florida è soltanto 975, Gore avrebbe vinto le elezioni presidenziali.