Gli esercizi sono riportati di seguito, prima elencando le domande, poi (a fine documento), gli svolgimenti. I numeri delle domande della prima sezione corrispondeono ai numeri delle risposte nella seconda.
Si hanno due campioni da 9 e 12 elementi. Il primo ha i valori 12, 11.7, 9.7, 11, 12.5, 12.7, 9.9, 11.4, 9.2; il secondo 12.5, 15.1, 14.4, 12.5, 14.9, 14.4, 13.2, 11.2, 12.6, 14.8, 13.9, 11.9. Assumendo che i campioni abbiano la stessa varianza, effettuare un test di Student a due lati senza utilizzare la funzione t.test().
Consideriamo il secondo campione dell’esercizio precedente. Risulta \(\bar y_2 = 13.45\). Verificare con la funzione t.test() l’ipotesi che il valore atteso della popolazione d’origine sia diverso da \(\mu_2 = 12.7\). Dato che nella fattispecie risulta \(\bar y_2 > \mu_2\), effettuare anche un test a un lato e discutere il risultato in termini di potenza del test.
Qual è il valore minimo di \(\mu_2\) che ci porta a rigettare \(H_0\) con una probabilità d’errore minore del 1%?
Considerare il dataset ggplot2::diamonds che riporta il prezzo di vendita di diamanti in funzione delle loro caratteristiche morfologiche.
Raggruppare la dimensione in carati (carat) in cinque gruppi di caratura. Costruire una tabella che riporti, per ciascun gruppo e per ciascuna categoria di taglio (cut), il minimo, il massimo e i tre quartili del prezzo di vendita, in colonne separate.
La tabella dati prod_defects.csv sul sito http://paolobosetti.quarto.pub/data.html#list-of-csv-files contiene dei dati di difettosità (come numero difetti per settimana) su 50 settimane di un impianto di produzione con tre macchine, M1, M2 e M3. Per le tre macchine i limiti massimi di difettosità sono, rispettivamente, 27, 32 e 100.
Creare una tabella che riporti, per ognuna delle 50 settimane, la difettosità in percentuale sul massimo per ciascuna macchina e realizzare il corrispondente grafico.
Per ogni macchina, aggiungere anche il calcolo della media progressiva, cioè il valor medio dei campioni dal primo all’i-esimo.
Ricordare che tutte le tabelle presenti in http://paolobosetti.quarto.pub/data.html possono essere caricate direttamente dalle funzioni read_table() o read_csv() passando l’URL generato come adas.utils::examples_url("nome_del_file.csv").
c1 <- c(12, 11.7, 9.7, 11, 12.5, 12.7, 9.9, 11.4, 9.2)
c2 <- c(12.5, 15.1, 14.4, 12.5, 14.9, 14.4, 13.2, 11.2, 12.6, 14.8, 13.9, 11.9)
c1.m <- mean(c1)
c2.m <- mean(c2)
n1 <- length(c1)
n2 <- length(c2)
Sp <- sqrt(((n1-1)*var(c1) + (n2-1)*var(c2)) / (n1 + n2 -2))
t0 <- abs(c1.m - c2.m) / (Sp * sqrt(1/n1 + 1/n2))
p.value <- pt(t0, n1 + n2 -2, lower.tail=FALSE) * 2
cat(paste("t0:", t0, "\n"))t0: 4.09946306243968 cat(paste("p-value:", p.value))p-value: 0.00061049441430277Per verifica:
t.test(c1, c2, var.equal = TRUE)
Two Sample t-test
data: c1 and c2
t = -4.0995, df = 19, p-value = 0.0006105
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.516249 -1.139306
sample estimates:
mean of x mean of y
11.12222 13.45000 A due lati:
t.test(c2, mu=12.7)
One Sample t-test
data: c2
t = 1.9921, df = 11, p-value = 0.07177
alternative hypothesis: true mean is not equal to 12.7
95 percent confidence interval:
12.62136 14.27864
sample estimates:
mean of x
13.45 A un lato:
(t <- t.test(c2, mu=12.7, alternative="greater", conf.level = 0.99))
One Sample t-test
data: c2
t = 1.9921, df = 11, p-value = 0.03588
alternative hypothesis: true mean is greater than 12.7
99 percent confidence interval:
12.42668 Inf
sample estimates:
mean of x
13.45 Adottando il test a un lato, il p-value risulta inferiore. Se ci fossimo dati una soglia per la probabilità di un falso positivo pari a \(\alpha=0.05\), in particolare, nel primo caso non avremmo potuto rigettare \(H_0\), mentre nel secondo sì.
Osservando l’intervallo di confidenza al 99%, valori di \(mu_2\) inferiori a 12.4266787 consentirebbero di rifiutare l’ipotesi nulla con una probabilità d’errore inferiore all’1%.
Provare a modificare il valore del parametro conf.level e osservare come l’unico risultato nell’output di t.test() è l’ampiezza dell’intervallo di confidenza.
diamonds %>%
mutate(size = cut(carat, breaks = 5), .after=carat) %>%
select(size:clarity, price) %>%
group_by(size,cut) %>%
summarise(
probs = seq(0, 100, length.out=5),
price = quantile(price),
.groups="keep"
) %>%
mutate(probs = paste0(probs, "%")) %>%
pivot_wider(names_from = probs, values_from = price) %>%
knitr::kable()| size | cut | 0% | 25% | 50% | 75% | 100% |
|---|---|---|---|---|---|---|
| (0.195,1.16] | Fair | 337 | 1715.00 | 2684.5 | 3785.25 | 10752 |
| (0.195,1.16] | Good | 327 | 896.50 | 2258.0 | 4009.50 | 17499 |
| (0.195,1.16] | Very Good | 336 | 773.00 | 1975.0 | 3937.00 | 18542 |
| (0.195,1.16] | Premium | 326 | 882.00 | 1826.5 | 4033.00 | 18279 |
| (0.195,1.16] | Ideal | 326 | 827.00 | 1389.0 | 2947.00 | 17590 |
| (1.16,2.12] | Fair | 2336 | 5334.00 | 7388.0 | 10913.00 | 18574 |
| (1.16,2.12] | Good | 3098 | 6686.50 | 8820.0 | 12236.00 | 18707 |
| (1.16,2.12] | Very Good | 2774 | 6945.25 | 9526.0 | 12762.75 | 18818 |
| (1.16,2.12] | Premium | 2360 | 6604.75 | 9384.5 | 12751.50 | 18795 |
| (1.16,2.12] | Ideal | 3168 | 7155.50 | 9761.0 | 12735.50 | 18806 |
| (2.12,3.09] | Fair | 5405 | 7069.00 | 12500.0 | 14975.00 | 18308 |
| (2.12,3.09] | Good | 9068 | 13751.50 | 14613.5 | 16894.25 | 18788 |
| (2.12,3.09] | Very Good | 6512 | 13782.00 | 15753.0 | 17161.00 | 18692 |
| (2.12,3.09] | Premium | 6357 | 13717.00 | 15857.5 | 17262.75 | 18823 |
| (2.12,3.09] | Ideal | 8709 | 14548.75 | 16427.0 | 17513.00 | 18791 |
| (3.09,4.05] | Fair | 9823 | 10745.50 | 11668.0 | 13816.00 | 15964 |
| (3.09,4.05] | Very Good | 15984 | 15984.00 | 15984.0 | 15984.00 | 15984 |
| (3.09,4.05] | Premium | 12300 | 15223.00 | 15223.0 | 16193.00 | 18701 |
| (3.09,4.05] | Ideal | 12545 | 12555.50 | 12566.0 | 12576.50 | 12587 |
| (4.05,5.01] | Fair | 17329 | 17673.50 | 18018.0 | 18274.50 | 18531 |
size:clarity: in dplyr è possibile selezionare più colonne dando un intervallo di nomi di colonne.summarise() può dare un warning suggerendo di usare reframe: durante il corso e durante l’esame NON utilizzate reframe(), perché essendo sperimentale non ne è garantita la disponibilità sui PC usati per l’esame.summarise() ritornano più di un valore, il sommario riporta altrettante righe. In questo caso come sommario utilizziamo quantile(), che ritorna 5 valori, e otteniamo quindi i valori desiderati, sebbene in un’unica colonna. Aggiungiamo la colonna perc per avere le corrispondenti percentuali.pivot_wider().Osserviamo il file di dati originale:
examples_url("prod_defects.csv") %>%
read_file() %>%
substr(1,250) %>% # solo i primi 250 caratteri
cat()# Difetti in produzione su 50 settimane
# S: numero settimana
# M1: conteggio difetti per macchina 1
# M2: conteggio difetti per macchina 2
# M3: conteggio difetti per macchina 3
S,M1,M2,M3
1,3,1,2
2,1,0,2
3,1,1,1
4,2,0,1
5,3,0,5
6,0,0,2
7,3,0,0
8,4,È evidente che il file contiene dei commenti inizianti con # e che i dati non sono tidy.
Importiamo quindi la tabella e trasformiamola in formato tidy:
prod <- examples_url("prod_defects.csv") %>%
read_csv(comment = "#") %>%
pivot_longer(-S, names_to="Macchina", values_to="Difetti") %>%
rename(Settimana=S)Visualizziamo i dati:
prod %>%
ggplot(aes(x=Settimana,y=Difetti, color=Macchina)) +
geom_point()
Creiamo una tabella di appoggio per i limiti:
limits <- tribble(
~Macchina, ~max,
"M1", 27,
"M2", 32,
"M3", 100
)Uniamo le due tabelle e calcoliamo il rapporto:
prod %>%
left_join(limits, by = join_by(Macchina)) %>%
mutate(`Dif.rel.`=(Difetti/max*100) %>% round(2)) %>%
ggplot(aes(x=Settimana, y=`Dif.rel.`, color=Macchina)) +
geom_line()
Per il calcolo della media progressiva possiamo usare una mappa, calcolando la media dei valori da 1 a \(i\):
examples_url("prod_defects.csv") %>%
read_csv(comment = "#") %>%
mutate(
M1.m = map_dbl(1:n(), \(i) mean(M1[1:i])),
M2.m = map_dbl(1:n(), \(i) mean(M2[1:i])),
M3.m = map_dbl(1:n(), \(i) mean(M3[1:i]))
) Rows: 50 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (4): S, M1, M2, M3
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# A tibble: 50 × 7
S M1 M2 M3 M1.m M2.m M3.m
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 3 1 2 3 1 2
2 2 1 0 2 2 0.5 2
3 3 1 1 1 1.67 0.667 1.67
4 4 2 0 1 1.75 0.5 1.5
5 5 3 0 5 2 0.4 2.2
6 6 0 0 2 1.67 0.333 2.17
7 7 3 0 0 1.86 0.286 1.86
8 8 4 0 0 2.12 0.25 1.62
9 9 2 0 1 2.11 0.222 1.56
10 10 2 1 3 2.1 0.3 1.7
# ℹ 40 more rowsOperare per colonne, però, non è molto pulito: supponiamo che in un futuro l’impianto cresca a 20 macchine, sarebbe necessario ripetere la stessa operazione più volte. È meglio risolvere il problema in formato tidy,
Per farlo partiamo dalla tabella tidy (prod). Non possiamo applicare direttamente una mappa, perché otterremo la media progressiva di tutti i difetti mescolando le varie macchine.
Possiamo però sfruttare group_by, che consente di operare su una tabella per gruppi:
prod %>%
group_by(Macchina) %>%
mutate(
mean = map_dbl(1:n(), ~ mean(Difetti[1:.]))
) %>%
slice_head(n=3)# A tibble: 9 × 4
# Groups: Macchina [3]
Settimana Macchina Difetti mean
<dbl> <chr> <dbl> <dbl>
1 1 M1 3 3
2 2 M1 1 2
3 3 M1 1 1.67
4 1 M2 1 1
5 2 M2 0 0.5
6 3 M2 1 0.667
7 1 M3 2 2
8 2 M3 2 2
9 3 M3 1 1.67 Nonostante slice_head(n=3) otteniamo 9 righe: come mai? è dovuto al fatto che la tabella è ancora raggruppata quindi ogni operazione si applica un gruppo alla volta, ottenendo tre righe per ognuno dei tre gruppi.
Il raggruppamento può essere disabilitato con ungroup().
Ora mettiamo tutto assieme, con l’attenzione di calcolare la media progressiva per la difettosità relativa (e non per il conteggio):
prod %>%
# Calcolo difettosità relativa
left_join(limits, by = join_by(Macchina)) %>%
mutate(`Dif.rel.`=(Difetti/max*100) %>% round(2)) %>%
# Calcolo le medie progressive per gruppi
group_by(Macchina) %>%
mutate(
mean = map_dbl(1:n(), ~ mean(`Dif.rel.`[1:.]))
) %>%
ungroup() %>%
# Grafico
ggplot(aes(x=Settimana, color=Macchina)) +
geom_line(aes(y=`Dif.rel.`)) +
geom_line(aes(y=mean), linetype=2)