On data science and IoT

Since several days have been reflecting on the deep connection between data science (or what we refer to predictive modeling) and IoT. IoT is commonly defined as (Rouse, 2019):

a network of interconnected computing devices, mechanical actuators and sensors able to exchange data between themselves without the need of human interaction.

It’s clear to me that the connection between this new network of things and data science is striking; In fact, I strongly believe the real revolution will come when these two branches of technology will finally be recognized as deeply related. I imagine a future where the data collected from the sensors will be transformed into insight and information by the machine learning algorithms and will automatically trigger a response in the physical world thanks to the physical actuators always connected on the Internet. 

Until now, data science has mostly focused on social network-generated data or Internet generated data (e.g. pictures, text mining on Twitter, etc); the insights that can be gathered from this kind of data is indeed limited in scope because no physical reaction can be triggered; or better, no improvement in efficiency can be triggered by using such data. On the contrary the data generated by the IoT world will pertain to the physical realm: think for example at the footfall in the city or the numbers of/the type of nutrient required by a crop field. All this data will be transmitted automatically and instantaneously over the Internet to algorithms able to predict and decide what to do based. This in turn will trigger a mechanical or chemical action inducing a response that is predetermined by humans using Machine Learning.

It is clear to me that the connection of the two technologies will be very important for humanity at large and it will be a multiplier of human capabilities in almost all fields of the physical realm.

Business Analytics con retaileR (parte II)

  1. Ecco la seconda parte dell’articolo sul nuovo pacchetto R retaileR. La prima parte potete trovarla qui.

Le Funzioni

closure.opp.cost

closure.opp.cost serve per stimare il costo di opportunità in funzione del tempo di chiusura dei locali di retail. Ossia, una stima del valore monetario perso a seconda dell’orario di chiusura.

I suoi argomenti sono:

Utilizzando il dataset preformattato, ecco un esempio:

closure_opp_cost("16:15:00", X= sales_august)

e l’output:

$Mean
[1] 10.35812

$Median
[1] 10.35

$StDev
[1] 8.392746

$`Late Sales`
[1] 10.30 18.90 6.30 26.70 14.60 10.60 0.05 24.18 10.40 9.60 14.10 14.90 5.10 0.00 0.00 0.00

In questo caso, chiudere il negozio alle 16:15:00 comporterebbe una perdita media di 10.35 €, con una varianza di 8.39 €. Inoltre, dal vettore Late.sales si può vedere chiaramente come 3 giorni non ci siano state vendite dopo le 16:15:00.

 segment_prod_line

Questa funzione è stata creata per l’analisi delle diverse famiglie di prodotti secondo un determinato segmento temporale.

Ossia, se la variabile temporale del dataset input è giornaliera, segment_prod_line produrrà le vendite medie giornaliere per la product line di interesse.

Vediamone il funzionamento:

specifichiamo la linea di prodotti come un vettore di stringhe:

prods <- c("Americano", "Espresso")

e successivamente evochiamo segment_prod_line

segment_prod_line(sales_august, prods)

ecco l’output:

[1] 31.56065

ciò significa che nei giorni che il dataset copre, abbiamo venduto in media 31.5 Americani ed Espressi al giorno.

Mean.items.times

Per ultima, una funzione di analisi temporale delle vendite. mean.items.times calcola le vendite medie (in volume) per ogni ora. Ossia crea un istogramma con l’unità temporale sull’asse delle ascisse e le vendite medie sull’asse delle ordinate.

Ecco il funzionamento (che ricordo, è ottimizzato per un dataset iZettle):

 mean.items.times(sales_august)

e l’output corrispondente:

       Time  Quantity
1  05:00:00  1.000000
2  06:00:00  9.000000
3  07:00:00 23.100000
4  08:00:00 30.695652
5  09:00:00 24.826087
6  10:00:00 17.695652
7  11:00:00 16.000000
8  12:00:00 11.173913
9  13:00:00 10.608696
10 14:00:00  8.590909
11 15:00:00  7.714286
12 16:00:00  6.000000
13 17:00:00  2.000000
14 18:00:00  1.000000

Business analytics con retaileR (parte I)

Perdonate la lunga assenza: sono stato molto impegnato in questi ultimi mesi. Tuttavia, per farmi perdonare, ho sviluppato una chicca tutta per voi:

un nuovissimo pacchetto R.

L’idea di creare un pacchetto R mi ha sempre affascinato. In fondo, è proprio grazie ai contributi Open Source che questo linguaggio continua ad arricchirsi. Tuttavia, mi è sempre mancato un ambito di applicazione, ossia non ho mai immaginato per cosa costruire il pacchetto. Per fortuna, negli ultimi mesi ho partecipato ad un progetto che mi ha introdotto al mondo del retail (il commercio Business to Consumer per intenderci) ed ai suoi peculiari problemi: da qui è nata l’idea di costruire una serie di funzioni che aiutino l’analisi dei dati in queso settore.

Il pacchetto, sviluppato in R, è chiamato retaileR e potete trovarlo sul mio Github. Nei prossimi mesi cercherò di renderlo disponibile sul CRAN, la repository ufficiale di tutti i pacchetti R.

Che cosa può fare retaileR? Retailer è nato con l’idea di creare un gruppo di funzioni utili per l’analisi dei dati di vendita. In questo modo, il processo decisionale verrebbe reso più semplice e con stime più efficaci. retaileR può aiutare nel decidere l’orario di chiusura di una filiale, organizzare lo staffing individuando gli orari di massimo/minimo volume e molte altre funzioni simili.

In questa serie di articoli illustrerò le funzioni e provvederò degli esempi su come utilizzare questo nuovo pacchetto.


I Dati

Il dataset di base viene considerato quello creato in automatico da iZettle, un sistema POS per il commercio al dettaglio. Qui sotto le prime righe di un dataset illustrativo (che userò negli esempi successivi) e che potrete trovare sempre in GitHub.

Carichiamo il dataset delle vendite fittizie relative ad un mese e visualizziamone le prime righe.

sales_august <- read.csv(“foo/sales_august.csv”)
head(sales_august)

Ecco l’output:

        Date                Time Receipt.number       Name    Variant Unit Quantity Price..GBP. Discount..GBP. Final.price..GBP.
1 2017-08-01 2017-08-20 06:19:00           7359  Americano              NA        1         2.5           0.00              2.50
2 2017-08-01 2017-08-20 06:19:00           7359   Espresso              NA        1         2.0           0.00              2.00
3 2017-08-01 2017-08-20 06:28:00           7360      Mocha              NA        1         3.1           1.55              1.55
4 2017-08-01 2017-08-20 06:43:00           7361 Flat White              NA        2         5.4           0.00              5.40
5 2017-08-01 2017-08-20 06:52:00           7362  Macchiato              NA        1         2.2           0.00              2.20
6 2017-08-01 2017-08-20 06:52:00           7362 Extra Shot Extra shot   NA        1         0.2           0.00              0.20

Le Funzioni

sales.format

questa funzione si propone come l’inizio di ogni analisi con retaileR poiché la sua funzione è formattare correttamente i dati da analizzare. In particolare, essa agisce sulle variabili di tempo (i.e. data e ora) formattandole in maniera uniforme.

Ecco un esempio:

sales_august <- sales.format(sales_sept)
[\code]

I suoi argomenti sono: X, date_var= “Date”, time_var= “Time”, format_time = ‘%H:%M’, dove:

  1. X: indica il data frame contenente i dati da analizzare.
  2. date_var è il nome, come stringa, della variabile contenente la data. Default a “Date”.
  3. time_var idem come date_var ma per il tempo (ossia l’ora). Default a “Time”.
  4. format_time come deve essere formattato il tempo. Per usare retaileR è necessario il formato ‘ora : minuto’ indicato appunto dal default ‘%H:%M’.
as.sales

l’analisi dei dati di vendita crea particolari problemi a livello della variabile temporale. Infatti, essendo una serie temporale, le vendite possono essere raggruppate seguendo tempi diversi: ore, giorni, settimane etc.

Inoltre, non è semplice capire se il di crescita sia positivo o negativo: la stagionalità e i movimenti casuali rendono molto spesso difficile questo compito.

Per evitare questi due problemi, la funzione as.sales crea una nuova classe di oggetti (in S4) che include una semplice regressione lineare e diversi raggruppamenti temporali.

i suoi argomenti sono:

  1. X: indica il data frame contenente i dati da analizzare.
  2. date_var è il nome, come stringa, della variabile contenente la data. Default a “Date”.
  3. time_var idem come date_var ma per il tempo (ossia l’ora). Default a “Time”.
  4. sales_var è il nome, come stringa, della variabile contenete il valore monetario di ogni singola vendita. Defaulta a “Final.price..GBP.”.
  5. format_time come deve essere formattato il tempo. Per usare retaileR è necessario il formato ‘ora : minuto’ indicato appunto dal default ‘%H:%M’.
  6. discount_var, nome colonna contenente gli sconti. Il default è “Discount..GBP.”.

Come output, la funzione produce un oggetto S4 sales, la cui struttura è la seguente:


ecco l'output:

 with 5 slots
  ..@ Daily           :'data.frame':	16 obs. of  2 variables:
  .. ..$ Date             : Date[1:16], format: "2017-08-01" "2017-08-02" "2017-08-03" "2017-08-04" ...
  .. ..$ Final.price..GBP.: num [1:16] 287 268 280 342 337 ...
  ..@ Weekly          :'data.frame':	3 obs. of  2 variables:
  .. ..$ Week             : num [1:3] 31 32 33
  .. ..$ Final.price..GBP.: num [1:3] 1514 2757 1565
  ..@ Monthly         :'data.frame':	1 obs. of  2 variables:
  .. ..$ Month            : num 8
  .. ..$ Final.price..GBP.: num 5836
  ..@ Discount        : num 39.8
  ..@ sales_funct_time:List of 12
  .. ..$ coefficients : Named num [1:2] -201711.1 11.6
  .. .. ..- attr(*, "names")= chr [1:2] "(Intercept)" "Date"
  .. ..$ residuals    : Named num [1:16] 9.79 -21.28 -20.64 29.43 13.04 ...
  .. .. ..- attr(*, "names")= chr [1:16] "1" "2" "3" "4" ...
  .. ..$ effects      : Named num [1:16] -1459 214.3 -18.6 30.4 13.1 ...
  .. .. ..- attr(*, "names")= chr [1:16] "(Intercept)" "Date" "" "" ...
  .. ..$ rank         : int 2
  .. ..$ fitted.values: Named num [1:16] 278 289 301 312 324 ...
  .. .. ..- attr(*, "names")= chr [1:16] "1" "2" "3" "4" ...
  .. ..$ assign       : int [1:2] 0 1
  .. ..$ qr           :List of 5
  .. .. ..$ qr   : num [1:16, 1:2] -4 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:16] "1" "2" "3" "4" ...
  .. .. .. .. ..$ : chr [1:2] "(Intercept)" "Date"
  .. .. .. ..- attr(*, "assign")= int [1:2] 0 1
  .. .. ..$ qraux: num [1:2] 1.25 1.27
  .. .. ..$ pivot: int [1:2] 1 2
  .. .. ..$ tol  : num 1e-07
  .. .. ..$ rank : int 2
  .. .. ..- attr(*, "class")= chr "qr"
  .. ..$ df.residual  : int 14
  .. ..$ xlevels      : Named list()
  .. ..$ call         : language lm(formula = as.formula(paste(sales_var, "~", date_var)), data = day)
  .. ..$ terms        :Classes 'terms', 'formula'  language Final.price..GBP. ~ Date
  .. .. .. ..- attr(*, "variables")= language list(Final.price..GBP., Date)
  .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:2] "Final.price..GBP." "Date"
  .. .. .. .. .. ..$ : chr "Date"
  .. .. .. ..- attr(*, "term.labels")= chr "Date"
  .. .. .. ..- attr(*, "order")= int 1
  .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. ..- attr(*, "response")= int 1
  .. .. .. ..- attr(*, ".Environment")=<environment: 0x102165000> 
  .. .. .. ..- attr(*, "predvars")= language list(Final.price..GBP., Date)
  .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "other"
  .. .. .. .. ..- attr(*, "names")= chr [1:2] "Final.price..GBP." "Date"
  .. ..$ model        :'data.frame':	16 obs. of  2 variables:
  .. .. ..$ Final.price..GBP.: num [1:16] 287 268 280 342 337 ...
  .. .. ..$ Date             : Date[1:16], format: "2017-08-01" "2017-08-02" "2017-08-03" "2017-08-04" ...
  .. .. ..- attr(*, "terms")=Classes 'terms', 'formula'  language Final.price..GBP. ~ Date
  .. .. .. .. ..- attr(*, "variables")= language list(Final.price..GBP., Date)
  .. .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. ..$ : chr [1:2] "Final.price..GBP." "Date"
  .. .. .. .. .. .. ..$ : chr "Date"
  .. .. .. .. ..- attr(*, "term.labels")= chr "Date"
  .. .. .. .. ..- attr(*, "order")= int 1
  .. .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. .. ..- attr(*, "response")= int 1
  .. .. .. .. ..- attr(*, ".Environment")=<environment: 0x102165000> 
  .. .. .. .. ..- attr(*, "predvars")= language list(Final.price..GBP., Date)
  .. .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "other"
  .. .. .. .. .. ..- attr(*, "names")= chr [1:2] "Final.price..GBP." "Date"
  .. ..- attr(*, "class")= chr "lm"

Come si può osservare, l’oggetto creato da as.sales contiene 5 sub-oggetti:

  •   3 raggruppamenti di tempo: daily, weekly, monthly.
  • Un valore dei discount totali nel dataset.
  • Una oggetto lm di regressione lineare, utile per comprendere il trend.

Breve digressione Natalizia sul confronto distribuzionale.

Recentemente al lavoro, ho dovuto scrivere un software che introduca della variabilità nei dati del Censo Scozzese. La motivazione di questa operazione è di proteggere la privacy dei cittadini, rendendo impossibile risalire agli utenti partendo dai dati.

Ovviamente, tale operazione, se fatta in maniera erronea, può portate a dei bias nei dati e inquinare tutte le potenziali analisi.

Dunque, come controllare che i dati processati siano simili ai dati originali?

In gergo statistico, fare questo procedimento significa controllare che le due serie di dati provengano dalla stessa distribuzione. Il test per eccellenza da utilizzare in questo caso è l’altisonante Kolmogorov-Smirnov. In generale, le comparazioni distribuzioni sono utilissime: ad esempio per controllare che misurazioni provenienti da due processi differenti possano essere analizzate usando gli stessi metodi. Insomma, per ipotizzare una somiglianza da due processi apparentemente diversi.

In R, il test di Kolmogorov-Smirnov può essere eseguito con la seguente funzione:

ks.test() {stats} 

 

Le assunzioni del modello lineare

In un articolo precedente abbiamo descritto la meccanica di un modello lineare. Come accennato, questi modelli sono costruiti su quattro assunzioni (because there’s no free lunch):

  • Normalità degli errori \epsilon_{i}
  • Deviazione standard degli errori costante
  • Indipendeza degli errori
  • Linearità della relazione tra la risposta e le variabili indipendenti

La violazione di questi requisiti causa conseguenze più o meno gravi sull’affidabilità del modello. Tuttavia, grazie a vari test, è possibile identificare e controllare se tali violazioni esistano. Nei paragrafi seguenti si illustrerà come eseguire questi test con R e la relativa matematica sottostante.

Da tenere a mente: gli errori di un modello sono definiti come:

\epsilon_{i} = y_{i} - \hat{y}_{i}

dove:

  • y_{i}: è il valore della risposta osservato.
  • \hat{y}_{i}: è il valore della risposta predetto dal modello.

Solitamente, per controllare queste assunzioni, gli errori devono essere standardizzati, ossia ridotti a seguire una distribuzione Normale con media 0 e dev standard uguale ad 1.

r_{i}^{s} = \frac{r_{i}}{\sqrt{s^{2}(1-h_{i,i})}}

h_{i,i} è la diagonale di una matrice, chiamata design matrix. Lasciamola da parte: l’importante è ricordarsi che gli errori standardizzati dovrebbero (se tutte le assunzioni sono rispettate) seguire una N(0,1).

Normalità

Gli errori dovrebbero essere approssimativamente distribuiti secondo una Normale con media zero e deviazione standard \sigma. Per controllare questa assunzione possiamo sia fare un test visivo che numerico.

Il test visivo consiste nel creare un QQ Plot, ossia raffigurare gli errori standardizzati contro un campione casuale di una N(0,1): se gli errori sono Normalmente distribuiti, essi devono posizionarsi sulla diagonale del piano:

Sinistra: errori distribuiti secondo N(0,1). Destra: normalità pesantemente violata.

Come test numerico invece, esistono numerose alternative (e siete liberi di procedere come meglio preferite). Per quanto mi riguarda uso (quasi) sempre il Shapiro-Wilk dove H0 è ‘Normalità rispettata’. In R:

shapiro.test(residuals(linearModel))

dove shapiro.test() è la funzione del test vera e propria e residuals(linearModel)) calcola gli errori per il nostro modello lineare.

La normalità tuttavia non è un’assunzione critica; violazioni minime possono essere comunque considerate accettabili.

Linearità

Per un modello lineare semplice, deve esistere una relazione lineare tra ogni variabile indipendente e quella dipendente. Per controllarla, basta un grafico della variabile in questione contro la dipendente: i punti devono posizionarsi in lungo una linea retta.

Omoschedasticità

Ossia varianza degli errori costante (chiamarla omoschedasticità però vi farà sembrare molto più professionali). Per controllarla, è necessario il test di Breusch-Pagan, dove H0 è ‘omoscedasticità rispettata’.

In R: ncvTest(linearModel)

Varianza degli errori costante è un’assunzione critica per un modello statistico: se non è rispettata può portarci a concludere che alcune variabili indipendenti siano utili per spiegare la dipendente, quando in realtà non lo sono. Ossia, ci può portare a trarre conclusioni errate sulle relazioni tra le variabili. Per risolvere questo problema, bisognerà abbandonare il modello lineare semplice per avventurarsi con qualche altro strumento (un GLS per esempio).

Indipendenza

Per indipendenza si intende che l’errore non deve essere correlato all’errore i + 1. Per controllarla, anche in questo caso, è possibile eseguire un test: il Durbin-Watson.

In R si esegue cosí durbinWatsonTest(linearModel). 

L’indipendenza è un’altra assunzione critica. Come per l’omoschedasticità, violazioni possono portare a erronee stime delle deviazioni standard, con conseguenti errori nei test di significatività. Un GLS potrebbe di nuovo essere una soluzione.

Regressione Lineare e diventare milionari

Regressione Lineare: predire le vendite

Ogni affarista che si rispetti dovrà prima o poi fare delle proiezioni di vendite: vuoi per ordinare il volume corretto dai fornitori, vuoi per dimostrare la validità del business ai finanziatori oppure per fare colpo con quanti affari farà il prossimo anno; le motivazioni sono molteplici e più o meno mondane. Ed ecco che il Data Scientist arriva in salvataggio! (E tornerà a casa con una bella parcella…)

Questo primo tutorial insegnerà proprio come, partendo da dati pre-esistenti, sia possible predire le future vendite di un business. Naturalmente, da buon accolito della statistica, non potrò esimermi da presentare la matematica sottostante: se ti annoia, salta alla sezione successiva dove mostrerò il lato pratico usando R.

1. Il Dataset

Il dataset che useremo si trova pre installato in R, ed è chiamato freeny. Di dimensioni 39×5 [39 righe/osservazioni x 5 colonne/variabili], rappresenta le vendite trimestrali di un non ben identificato business. Le variabili presenti sono:

  1. y: le vendite. Notate la notazione: in statistica/Data Science,  identifica sempre la variabile da predire, anche chiamata variabile dipendente risposta (in inglese, dependent variable/response).
  2. lag.quarterly.revenue: le vendite del trimestre precedente.
  3. price.index: livello dei prezzi dei beni venduti.
  4. income.level: livello degli stipendi della clientela.
  5. market.potential: indice dell’attrattività del mercato.

Ecco una vista di come appare:

Il dataset freeny.

Nota: per creare modelli di forecasting utilizzeremo sempre dataset in forma matriciale.

2. La Regressione Lineare (LM)

la matematica no!

Come accennato precedentemente, il nostro obbiettivo è quello di predire le vendite per un trimestre futuro. Qual’è il ragionamento da seguire? Noi, esseri mortali, non conosciamo le vendite future dato che non si sono ancora verificate. Tuttavia, abbiamo a nostra disposizione una serie di altre variabili di cui conosciamo perfettamente il valore, anche futuro. In questo caso ad esempio, conosciamo le vendite passate (lag.quarterly.revenue), possiamo decidere che prezzo mettere (price.index), sappiamo il livello degli stipendi (income.level) tramite l’Eurostat e conosciamo anche il potenziale di un certo mercato (market.potential), magari tramite un sondaggio. Qual’è dunque il procedimento più logico per raggiungere il nostro obbiettivo? Spiegare le vendite come una funzione delle variabili a noi note; in altre parole, cercare una relazione tra le vendite e tutto il resto.

In Statistica, questo tipo di relazione, viene chiamato regressione. La forma funzionale più semplice è la seguente:

y_{i} = \beta_{0} + \beta_{1}x_{i} + \epsilon_{i}

Dove:

  • y_{i} : è il valore per le vendite (variabile dipendente) all’osservazione i, ossia alla riga i
  • x_{i} : è il valore per una variabile indipendente (una di quelle che usiamo per predire) all’osservazione i
  • \beta_{0} : è un parametro, anche chiamato intercetta (vedremo tra poco il motivo)
  • \beta_{1} :  è un altro parametro, anche chiamato pendenza, che rappresenta il cambio di per ogni unità di differenza in x. Ossia se \beta_{1} = 2, e aumenta di 1, allora aumenterà di 2.
  • \epsilon_{i} : i residuali, la componente casuale, che riamane insperata dal modello.

Gli amanti della geometria noteranno che tale funzione non è altro che una retta (da qui i termini ‘intercetta’ e ‘pendenza’). E infatti, quello che cerchiamo di fare è proprio questo: cercare di creare una retta che si adatti il meglio possibile ai dati.

Ad esempio:

Vendite per livello stipendi

visualizza la relazione esistente tra le vendite e il livello degli stipendi. Come possiamo vedere, essa non è una retta perfetta: questo perché c’è ovviamente una componente aleatoria nella relazione. Costruendo una retta perfetta che si adattati al meglio a questi dati, andremo proprio a separare la relazione ‘vera’ dalla componete aleatoria; separeremo il segnale dal rumore di fondo, che è casuale. Ecco il risultato di questo procedimento:

Linea di regressione (in rosso)

Per costruire questa retta dobbiamo ‘aggiustare’ i due parametri \beta_{0}, \beta_{1}. In che modo? Minimizzando la distanza tra ogni punto e la retta!

Ossia usando il metodo dei least squaresi quadrati minimi; minimizzando questo tizio qui:

\sum\limits_{i=1}^{n}(y_{i} - (\beta_{0} + \beta{1}x_{i}))^{2}

Questo caso con una sola variabile indipendente può essere poi generalizzato a più variabili:

y_{i} = \beta_{0} + \beta_{1}x_{1i} + \beta_{2}x_{2i} + \cdots + \beta_{p}x_{pi} + \epsilon_{i}

La logica è sempre la stessa, anche se ora è impossibile visualizzare la relazione con un grafico (il massimo possibile sono 3 dimensioni. Per ovvie ragioni).

Naturalmente, come ogni elemento statistico, anche i parametri della regressione avranno una componente di variabilità, rappresentata dalla loro varianza e deviazione standard.

Inoltre, i modelli di regressione lineare semplice sono basati su quattro assunzioni:

  • Normalità degli errori \epsilon_{i}.
  • Deviazione standard degli errori \epsilon_{i} costante.
  • Linearità della relazione tra le variabili.
  • Indipendenza degli errori.

Data la vastità dell’argomento, delle assunzioni parleremo più nel dettaglio in un articolo successivo.

3. In R

Finalmente la pratica.

La regressione lineare semplice, in R, richiede ben pochi sforzi. La funzione da utilizzare è lm() nel pacchetto di {stats} (che è già preinstallato).

Per creare il modello basta dunque scrivere :

attach(freeny) # Per caricare il dataset
View(freeny) # Per visionare il dataset in forma matriciale
modello = lm(y~., data = freeny) # Per creare il modello

Il modello cosí creato avrà come variabile dipendente e tutte le altre (il ‘.’) come indipendenti.

Per visionare il contenuto, ovvero i parametri del modello, la loro variabilità e i relativi p-value, bisogna utilizzare la funzione summary():

summary(modello)

Ecco l’output:

Output del LM

Come possiamo vedere, la variabile log.quarterly.revenue, quando tutte le altre variabili sono presenti, non risulta avere una relazione statisticamente significativa con le vendite.