Dati sui succhi di frutta

Consideriamo dei dati sulle preferenze di consumatori di due marche (MM e CH) di succhi di frutta derivanti da \(n = 1070\) acquisti che includono il succo di frutta. La previsione della scelta del consumatore è fatta sulla base delle variabili seguenti: i prezzi dei due prodotti priceCH e priceMM, gli sconti applicati discountCH and discountMM, un indicatore della fedelta per il prodotto MM loyaltyMM, un indicatore sella settimana di acquisto week, e il supermercato dove è stato fatto l’acquisto, store. L’analisi seguente è una traduzione quasi letterale di Azzalini e Scarpa (2012). Data analysis and data mining. Wiley, New York.

Variable     Description 

choice       pre-chosen brand (factor, with 2 levels)
id.cust      customer identification 
week         identifier of week of purchase 
priceCH      reference price for  brand CH (USD)
priceMM      reference price for  brand MM (USD)
discountCH   discount applied to  product CH (USD)
discountMM   discount applied to product MM (USD)
loyaltyCH    loyalty indicator for  product CH 
loyaltyMM    loyalty indicator for  product MM 
store        store identifier (factor, with 5 levels)

NOTA: La variabile loyaltyMM è costruita partendo dal valore 0.5 e aumentandolo per ogni acquisto dello stesso consumatore, con un valore pari al 20% della differenza tra il valore corrente e 1 se l’acquirente sceglie MM e diminuendolo del 20% della differenza tra il valore corrente e 0 se l’acquirente sceglie CH. La variabile loyaltyCH è dato da 1-loyaltyMM.

Lettura dei dati

rm(list = ls())
juice =  read.table("~/Dropbox/Public/gmm/static/R/juice.data.txt", head=TRUE) 
#juice =  read.table("http://labdisia.disia.unifi.it/gmm/downloads/files/juice.data.txt", head=TRUE) 
juice[,"store"] = factor(juice[,"store"])
v = c(1,3,6:9,13,21)
juice = juice[,v]

Selezioniamo una porzione aleatoria comprendente i 3/4 dei dati totali per adattare il modello. I dati restanti sono usati per valutare i risultati.

set.seed(123)
n = nrow(juice)
n1 = round(n*0.75)
n2= n-n1
permuta= sample(1:n,n)
training = sort(permuta[1:n1])
#
juice1 = juice[training,]
juice2 = juice[-training,]

Analisi preliminare

Box- delle variabili continue classificando per le due marche e un grafico a barre delle percentuali di preferenza dalla marca MM in ciascun supermercato.

par(mfrow=c(2,3))
K= length(v)
for(k in 2:(K-1))
boxplot(juice1[,k]~choice, col= "orange",   ylab=names(juice1)[k], data = juice1)

par(mfrow = c(1,1))
freq= with(juice1, table(choice, store))
freq= 100*freq/outer(rep(1,2),apply(freq,2,sum))
barplot(freq[2,], width=1, space=1.25, col="RoyalBlue", ylim=c(0,100),
ylab="Percentage of choosing MM", xlab="Shop")
abline(h=100* mean(juice1$choice=="MM") , lty=2, col="red")

Modello logistico

m1 = glm(choice ~ week + priceCH + priceMM + discountCH + discountMM +  loyaltyMM + store, data=juice1, family=binomial)
summary(m1)

Call:
glm(formula = choice ~ week + priceCH + priceMM + discountCH + 
    discountMM + loyaltyMM + store, family = binomial, data = juice1)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.749  -0.567  -0.250   0.555   2.779  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.81611    2.06284   -1.85  0.06432 .  
week        -0.00173    0.01309   -0.13  0.89486    
priceCH      4.43557    2.11851    2.09  0.03628 *  
priceMM     -3.70577    1.00809   -3.68  0.00024 ***
discountCH  -3.64867    1.14489   -3.19  0.00144 ** 
discountMM   2.09456    0.50162    4.18    3e-05 ***
loyaltyMM    5.86455    0.44978   13.04  < 2e-16 ***
store1       0.55120    0.31550    1.75  0.08062 .  
store2       0.65630    0.28592    2.30  0.02171 *  
store3       0.57354    0.36819    1.56  0.11930    
store4       0.03861    0.42014    0.09  0.92677    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1077.20  on 801  degrees of freedom
Residual deviance:  631.63  on 791  degrees of freedom
AIC: 653.6

Number of Fisher Scoring iterations: 5

Riduzione del modello

m2 = update(m1, .~ .-week, data=juice1)
summary(m2)

Call:
glm(formula = choice ~ priceCH + priceMM + discountCH + discountMM + 
    loyaltyMM + store, family = binomial, data = juice1)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.758  -0.566  -0.252   0.553   2.779  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -3.8127     2.0630   -1.85  0.06458 .  
priceCH       4.2411     1.5233    2.78  0.00537 ** 
priceMM      -3.7444     0.9652   -3.88  0.00010 ***
discountCH   -3.6952     1.0892   -3.39  0.00069 ***
discountMM    2.0818     0.4918    4.23  2.3e-05 ***
loyaltyMM     5.8683     0.4489   13.07  < 2e-16 ***
store1        0.5433     0.3095    1.76  0.07923 .  
store2        0.6513     0.2834    2.30  0.02155 *  
store3        0.5927     0.3386    1.75  0.08001 .  
store4        0.0548     0.4021    0.14  0.89168    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1077.20  on 801  degrees of freedom
Residual deviance:  631.64  on 792  degrees of freedom
AIC: 651.6

Number of Fisher Scoring iterations: 5

Test del rapporto di veosimiglianza (1 gl)

anova(m2, m1, test = "Chisq")
Analysis of Deviance Table

Model 1: choice ~ priceCH + priceMM + discountCH + discountMM + loyaltyMM + 
    store
Model 2: choice ~ week + priceCH + priceMM + discountCH + discountMM + 
    loyaltyMM + store
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       792        632                     
2       791        632  1   0.0175     0.89

Test con più gradi di libertà

m3 = update(m1, . ~ .- week - store, data=juice1)
summary(m3)

Call:
glm(formula = choice ~ priceCH + priceMM + discountCH + discountMM + 
    loyaltyMM, family = binomial, data = juice1)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-2.713  -0.584  -0.255   0.564   2.734  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -3.366      1.804   -1.87   0.0620 .  
priceCH        4.131      1.286    3.21   0.0013 ** 
priceMM       -3.731      0.912   -4.09  4.3e-05 ***
discountCH    -4.193      1.022   -4.10  4.1e-05 ***
discountMM     2.017      0.464    4.35  1.4e-05 ***
loyaltyMM      6.119      0.432   14.15  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1077.20  on 801  degrees of freedom
Residual deviance:  639.43  on 796  degrees of freedom
AIC: 651.4

Number of Fisher Scoring iterations: 5
anova(m3, m1, test = "Chisq")
Analysis of Deviance Table

Model 1: choice ~ priceCH + priceMM + discountCH + discountMM + loyaltyMM
Model 2: choice ~ week + priceCH + priceMM + discountCH + discountMM + 
    loyaltyMM + store
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       796        639                     
2       791        632  5      7.8     0.17

Tavole di errata classificazione

Ora applichiamo il modello m2 ai test data per classificare gli acquisti rimanenti e valutare la capacità predittiva del classificatore appena trovato e costruiamo la matrice di confusione:

\[ \begin{array}{cccc} Classificati &Effettivi& - & + & Tot \\ \hline - && n_{00} & n_{01} & n_{0\cdot}\\ + && n_{10} & n_{11} & n_{1\cdot}\\ \hline Tot && n_{\cdot 0} & n_{\cdot 1} & n\\ \end{array} \]

`matrice.confusione` = function(previsti, osservati){
  a <-  table(previsti, osservati)
  err.tot <- 1-sum(diag(a))/sum(a)
  fn <- a[1,2]/(a[1,2]+a[2,2])
  fp <- a[2,1]/(a[1,1]+a[2,1])
  print(a)
  cat("total error: ", format(err.tot),"\n")
  cat("false positives/negatives: ",format(c(fp, fn)),"\n")
  invisible(a)
}

p2 = predict(m2, newdata=juice2, type="response")
matrice.confusione(p2>0.5, juice2[,1])
        osservati
previsti  CH  MM
   FALSE 150  23
   TRUE   19  76
total error:  0.1567 
false positives/negatives:  0.1124 0.2323

I casi classificati correttamente e l’errore totale sono \[ (150 + 76)/268 = 0.843, \quad (19 + 23)/268 = 0.157 \] Se chiamiamo “positivi” quelli che scelgono MM,

  • i falsi positivi sono gli acquirenti classificati MM mentre erano acquirenti di CH.
  • i falsi negativi sono gli acquirenti classificati CH mentre erano acquirenti di MM.

Le rispettive probabilità sono stimate come
\[ FP = 19/(150+19) = 0.1124, \qquad FN = 23/(23+76) = 0.2323 \]

Curva ROC

Per valutare l’adeguatezza del classificatore si può usare anche la curva ROC (receiver operating characteristic). È stata introdotta durante la II guerra mondiale nell’ambito della teora delle comunicazioni (capacità di individuare i segnali radar) e applicata poi nel controllo di qualità e in statistica medica.

Le proporzioni di falsi positivi e falsi negativi varia al variare del threshold \(1/2\) usato per effettuare la classificazione. Se la regola è \[ \text{Classifica come MM se } \hat P(\text{MM}) > 1/2 \] si ottiene la tavola di confusione precedente. Se il threshold viene fatto variare la tavola cambia come pure le proporzioni di falsi positivi e negativi.

  • Si dice specificità la proporzione di classificati negativi sul totale di effettivi negativi

  • Si dice sensibilità la proporzione di classificati positivi sul totale di effettivi positivi

La curva ROC è una spezzata che collega i punti con coordinate \[ (1- \text{specificità}, \text{sensibilità}) \] cioè le proporzioni di falsi positivi e di veri positivi, al variare del threshold. Nota che sono l’equivalente del livello e della potenza del test.

g = as.numeric(juice2[,1]=="MM")
require(pROC)
plot(roc(g, p2), legacy.axes = TRUE)

Call:
roc.default(response = g, predictor = p2)

Data: p2 in 169 controls (g 0) < 99 cases (g 1).
Area under the curve: 0.914
# xlab = "Average false positive rate", ylab = "Average true positive rate"
abline(v = 1-0.1124, col = "red", lty = 2)
abline(h = 1- 0.2323, col = "red", lty = 2)

Predizione col classificatore m3

p3 = predict(m3, newdata=juice2, type="response")
matrice.confusione(p3>0.5, juice2[,1])
        osservati
previsti  CH  MM
   FALSE 147  25
   TRUE   22  74
total error:  0.1754 
false positives/negatives:  0.1302 0.2525 
plot(roc(g, p3), xlab = "Average false positive rate", ylab = "Average true positive rate")


Call:
roc.default(response = g, predictor = p3)

Data: p3 in 169 controls (g 0) < 99 cases (g 1).
Area under the curve: 0.913

Un modo per visualizzare gli effetti non-lineari

require(effects)
plot(allEffects(m3), main = NULL)

Esercizi 9

  1. Verificare che l’informzione attesa non è invariante rispetto alle riparametrizzazioni. Usare come esempio il modello di campionamento iid da una densità esponenziale \[ f(y; \lambda) = \theta e^{-\theta y}, \quad y > 0, \theta >0 \] e la riparametrizzazione \(\beta = 1/\theta\). Calcolate gli errori standard per \(\hat \theta\) e \(\hat \beta\).

  2. Si consideri il data frame contr.small. Il data frame si ottiene con il comando:

source("http://labdisia.disia.unifi.it/gmm/downloads/files/contr.small.R")
head(contr.small)
   use age indian urban ed
3    0  35      0     0  8
6    0  37      0     0 16
13   1  41      0     0  5
16   1  42      1     0  1
17   1  37      1     0  0
18   1  36      1     0  5

I dati riguardano lo studio della fecondità e la contraccezione (n. di figli rilevato durante la World Fertility Survey) nelle isole Fiji nel 1974. È stato considerato un campione di 963 donne fertili, sposate e non incinte di età compresa tra 35 e 44 anni e si sono rilevate le variabili

Adattare un modello logistico per studiare la dipendenza della probabilità di usare un contraccettivo in funzione delle variabili esplicative. Valutare l’effetto del livello di istruzione al netto delle altre variabili. Valutare se esistono effetti di tipo non lineare o interattivi.

Soluzione

  1. L’informazione attesa per \(\theta\) per l’esponenziale si trova così. \[ l_i(\theta) = log \theta - Y_i \theta, \quad U_i = \frac{1}{\theta} - Y_i \] Quindi \(U = \frac{n}{\theta} - \sum_i Y_i\). Dunque \[ I(\theta) = V(U) = \sum_i V(Y_i) = n/\theta^2. \]

L’informazione attesa per \(\beta\) analogamente si trova con \[ l_i(\beta) = -\log (beta) - Y_i /\beta , \qquad U_i = -1/\beta + Y_i/\beta^2 \] Perciò \[ V(\sum_i U_i) = + \frac{1}{\beta^4} V(\sum_i Y_i ) = \frac{n\beta^2}{\beta^4} = n/\beta^2 \]

In generale se \(\beta = 1/\theta\) e \(\theta = g(\beta) = 1/\beta\) \[ I(\beta) = (g'(\beta))^2 I(\theta) = \frac{1}{\beta^4}I(\theta) \]

  1. Istruzioni
ml = glm(use ~ age + ed + indian + urban, binomial, data = contr.small)
summary(ml)

Call:
glm(formula = use ~ age + ed + indian + urban, family = binomial, 
    data = contr.small)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.848  -1.108   0.764   0.943   1.436  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.5621     1.0017   -1.56    0.119    
age           0.0249     0.0250    1.00    0.319    
ed            0.0541     0.0249    2.17    0.030 *  
indian        1.1683     0.1662    7.03  2.1e-12 ***
urban         0.3676     0.1533    2.40    0.016 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1291.9  on 962  degrees of freedom
Residual deviance: 1220.2  on 958  degrees of freedom
AIC: 1230

Number of Fisher Scoring iterations: 4