Il set di dati ha 14 campi.
- Codice VCAA: ID amministrativo per ciascun codice
- Nome della scuola
- Settore – G = Governo, C = Cattolico e I = Indipendente
- Località o sobborgo
- Totale studenti che completano l’anno 12
- Conmittenti in pista
- Rispondenti
- Il resto delle colonne rappresenta la percentuale di intervistati per ciascun risultato
Ai fini del nostro studio trasversale siamo interessati alle proporzioni dei risultati per settore e quindi dobbiamo convertire questo ampio frame di dati in un formato lungo.
df_long <- df |>
mutate(across(5:14, ~ as.numeric(.x)), #convert all numeric fields captured as characters
across(8:14, ~ .x/100 * respondants), #calculate counts from proportions
across(8:14, ~ round(.x, 0)), #round to whole integers
respondants = bachelors + deferred + tafe + apprentice_trainee + employed + looking_for_work + other) |> #recalculate total respondents |>
filter(sector != 'A') |> #Low volume
select(sector, school_name, 7:14) |>
pivot_longer(c(-1, -2, -3), names_to = 'outcome', values_to = 'proportion') |>
mutate(proportion = proportion/respondants)
Analisi esplorativa dei dati
Visualizziamo e riassumiamo brevemente il nostro set di dati. Il settore governativo ne ha 157, l’Independent ne ha 57 e il Catholic ha 50 scuole elencate in questo set di dati.
df |>
mutate(sector = fct_infreq(sector)) |>
ggplot(aes(sector)) +
geom_bar(aes(fill = sector)) +
geom_label(aes(label = after_stat(count)), stat = 'count', vjust = 1.5) +
labs(x = 'Sector', y = 'Count', title = 'Count of Schools by Sector', fill = 'Sector') +
scale_fill_viridis_d(begin = 0.2, end = 0.8) +
theme_ggdist()
df_long |>
ggplot(aes(sector, proportion, fill = outcome)) +
geom_boxplot(alpha = 0.8) +
geom_point(position = position_dodge(width=0.75), alpha = 0.1, color = 'grey', aes(group = outcome)) +
labs(x = 'Sector', y = 'Proportion', fill = 'Outcome', title = 'Boxplot of Respondant Proportions by Sector and Outcome') +
scale_fill_viridis_d(begin = 0.2, end = 0.8) +
theme_ggdist()
Le distribuzioni raccontano una storia interessante. I risultati delle lauree mostrano la maggiore variabilità in tutti i settori, con le scuole indipendenti che hanno la percentuale mediana di studenti più elevata che sceglie di proseguire gli studi universitari. È interessante notare che le scuole statali hanno la percentuale media più alta di studenti che trovano lavoro dopo aver completato la scuola superiore. Tutti gli altri risultati non variano troppo: lo rivisiteremo presto.
Modellazione statistica: regressione della probabilità beta
Ci concentriamo sulle proporzioni degli studenti per scuola e sui loro risultati dopo il diploma di scuola superiore. Una probabilità beta è la nostra migliore scommessa in questi casi. brms è un pacchetto brillante di Buerkner per sviluppare modelli bayesiani. Il nostro obiettivo è modellare la distribuzione delle proporzioni per risultato e settore.
La regressione beta modella due parametri, μ e φ. μ rappresenta la proporzione media e φ è una misura di precisione (o varianza).
Il nostro primo modello è rappresentato di seguito, nota che stiamo già iniziando con un’interazione tra Settore e Risultato perché questa è la domanda a cui vogliamo che il modello risponda, e quindi questo è un modello ANOVA.
Chiediamo al modello di creare termini Beta individuali per ciascuna combinazione di Settore e Risultato, con un termine φ raggruppato, o di modellare medie proporzionali diverse con la stessa varianza. Ci aspettiamo che il 50% di queste proporzioni sia compreso tra (-3, 1) sulla scala logit o (0,01, 0,73) come proporzioni. Questo è sufficientemente ampio ma informato in anticipo. La stima Phi globale è un numero positivo, quindi utilizziamo una distribuzione gamma sufficientemente ampia.
# Modelling Proportion with Sector:Outcome Interaction term using Beta - Pooled Phi termm3a <- brm(
proportion ~ sector:outcome + 0,
family = Beta,
data = df_long |> mutate(proportion = proportion + 0.01), # Beta regression can't take outcomes that are 0 so we fudge here by adding tiny increment
prior = c(prior(normal(-1, 2), class = 'b'),
prior(gamma(6, 1), class = 'phi')),
iter = 8000, warmup = 1000, cores = 4, chains = 4, seed = 246, save_pars = save_pars(all = T),
control = list(adapt_delta = 0.99, max_treedepth = 15)
) |>
add_criterion(c('waic', 'loo'), moment_match = T)
summary(m3a)
Tieni presente che nella configurazione del modello abbiamo aggiunto un incremento dell’1% a tutte le proporzioni: questo perché la regressione Beta non gestisce valori zero. Abbiamo tentato di modellarlo con un beta gonfiato pari a zero, ma ci è voluto molto più tempo per convergere.
Allo stesso modo, possiamo modellare senza un phi aggregato, questo intuitivamente ha senso dato ciò che abbiamo osservato nelle distribuzioni di cui sopra, c’è variazione tra ciascuna combinazione di settore e risultato. Il modello è definito di seguito.
m3b <- brm(
bf(proportion ~ sector:outcome + 0,
phi ~ sector:outcome + 0),
family = Beta,
data = df_long |> mutate(proportion = proportion + 0.01),
prior = c(prior(normal(-1, 2), class = 'b')),
iter = 8000, warmup = 1000, cores = 4, chains = 4, seed = 246, save_pars = save_pars(all = T),
control = list(adapt_delta = 0.99)
) |>
add_criterion(c('waic', 'loo'), moment_match = T)summary(m3b)
Diagnostica e confronto dei modelli
Con due modelli in mano, confronteremo ora la loro accuratezza predittiva fuori campione utilizzando la stima bayesiana LOO della densità predittiva puntuale logaritmica prevista (elpd_loo).
comp <- loo_compare(m3a, m3b)print(comp, simplify = F)
In poche parole, maggiore è il valore di esclusione di un punto logaritmico previsto, maggiore è l’accuratezza predittiva sui dati invisibili. Questo ci dà un bene parente misura dell’accuratezza del modello tra modelli. Possiamo verificarlo ulteriormente completando un controllo predittivo a posteriori, un confronto tra valori osservati e simulati. Nel nostro caso, il modello m3b è il modello migliore dei dati osservati.
alt_df <- df_long |>
select(sector, outcome, proportion) |>
rename(value = proportion) |>
mutate(y = 'y',
draw = 99) |>
select(sector, outcome, draw, value, y)sim_df <- expand_grid(sector = c('C', 'I', 'G'),
outcome = unique(df_long$outcome)) |>
add_predicted_draws(m3b, ndraws = 1200) |>
rename(value = .prediction) |>
ungroup() |>
mutate(y = 'y_rep',
draw = rep(seq(from = 1, to = 50, length.out = 50), times = 504)) |>
select(sector, outcome, draw, value, y) |>
bind_rows(alt_df)
sim_df |>
ggplot(aes(value, group = draw)) +
geom_density(aes(color = y)) +
facet_grid(outcome ~ sector, scales = 'free_y') +
scale_color_manual(name = '',
values = c('y' = "darkblue",
'y_rep' = "grey")) +
theme_ggdist() +
labs(y = 'Density', x = 'y', title = 'Distribution of Observed and Replicated Proportions by Sector and Outcome')
Il modello m3b, data la sua varianza non raggruppata o termine phi, è in grado di catturare meglio la variazione nelle distribuzioni proporzionali per settore e risultato.
ANOVA: stile bayesiano
Ricordiamo che la nostra domanda di ricerca mira a cercare di capire se ci sono differenze nei risultati tra i settori e in che misura. Nelle statistiche frequentiste, potremmo usare ANOVA, un approccio con differenza di medie tra gruppi. Il punto debole di questo approccio è che i risultati forniscono una stima e un intervallo di confidenza, senza senso di incertezza su queste stime, e un valore p controintuitivo che indica se la differenza nelle medie è statisticamente significativa o meno. No grazie.
Di seguito generiamo una serie di valori attesi per ogni settore e interazione con la combinazione di risultati, quindi utilizziamo l’eccellente funzione tidybayes::compare_levels() che esegue il lavoro pesante. Calcola la differenza nelle medie a posteriori tra i settori per ciascun risultato. Abbiamo escluso l’esito “altro” per brevità.
new_df <- expand_grid(sector = c('I', 'G', 'C'),
outcome = c('apprentice_trainee', 'bachelors', 'deferred', 'employed', 'looking_for_work', 'tafe'))epred_draws(m3b, newdata = new_df) |>
compare_levels(.epred, by = sector, comparison = rlang::exprs(C - G, I - G, C - I)) |>
mutate(sector = fct_inorder(sector),
sector = fct_shift(sector, -1),
sector = fct_rev(sector)) |>
ggplot(aes(x = .epred, y = sector, fill = sector)) +
stat_halfeye() +
geom_vline(xintercept = 0, lty = 2) +
facet_wrap(~ outcome, scales = 'free_x') +
theme_ggdist() +
theme(legend.position = 'bottom') +
scale_fill_viridis_d(begin = 0.4, end = 0.8) +
labs(x = 'Proportional Difference', y = 'Sector', title = 'Differences in Posterior Means by Sector and Outcome', fill = 'Sector')
In alternativa possiamo riassumere tutte queste distribuzioni con una tabella ordinata per una più facile interpretazione e un intervallo credibile dell’89%.
marg_gt <- epred_draws(m3b, newdata = new_df) |>
compare_levels(.epred, by = sector, comparison = rlang::exprs(C - G, I - G, C - I)) |>
median_qi(.width = .89) |>
mutate(across(where(is.numeric), ~round(.x, 3))) |>
select(-c(.point, .interval, .width)) |>
arrange(outcome, sector) |>
rename(diff_in_means = .epred,
Q5.5 = .lower,
Q94.5 = .upper) |>
group_by(outcome) |>
gt() |>
tab_header(title = 'Sector Marginal Distributions by Outcome') |>
#tab_stubhead(label = 'Sector Comparison') |>
fmt_percent() |>
gtExtras::gt_theme_pff()
Ad esempio se dovessimo riassumere un confronto in una relazione formale potremmo scrivere quanto segue.
Gli studenti delle scuole statali hanno meno probabilità di iniziare i corsi di laurea rispetto ai loro omologhi delle scuole cattoliche e indipendenti dopo il completamento del VCE.
In media, il 42,5% (tra il 39,5% e il 45,6%) degli studenti delle scuole statali, il 53,2% (tra il 47,7% e il 58,4%) degli studenti delle scuole cattoliche e il 65% (tra il 60,1% e il 69,7%) degli studenti delle scuole indipendenti iniziano gli studi universitari dopo il completamento degli studenti del 12° anno.
Esiste una probabilità a posteriori dell’89% che la differenza tra l’iscrizione agli studenti universitari cattolici e quelli governativi sia compresa tra il 5,6% e il 15,7% con una media del 10,7%. Inoltre, la differenza tra l’iscrizione agli studenti universitari indipendenti e governativi è compresa tra il 17,8% e il 27% con una media del 22,5%.
Queste differenze sono sostanziali e c’è una probabilità del 100% che queste differenze non siano pari a zero.
Sommario e conclusione
In questo articolo abbiamo dimostrato come modellare dati proporzionali utilizzando una funzione di verosimiglianza Beta utilizzando il modello bayesiano e quindi eseguire l’ANOVA bayesiana per esplorare le differenze nei risultati proporzionali tra i settori.
Non abbiamo cercato di creare una comprensione causale di queste differenze. Si può immaginare che ci siano diversi fattori che influenzano il rendimento dei singoli studenti, il background socioeconomico, il livello di istruzione dei genitori, oltre all’impatto a livello scolastico, alle risorse, ecc. Si tratta di un’area incredibilmente complessa della politica pubblica che spesso può impantanarsi in zero- retorica sommaria.
Personalmente, sono la prima persona nella mia famiglia allargata a frequentare e completare l’istruzione terziaria. Ho frequentato una scuola superiore pubblica di fascia media e ho ottenuto punteggi ragionevolmente buoni per ottenere l’ingresso nella mia prima preferenza. I miei genitori mi hanno incoraggiato a proseguire gli studi, decidendo entrambi di lasciare la scuola all’età di 16 anni. Sebbene questo articolo fornisca la prova che la differenza tra scuole governative e non governative è reale, è di natura puramente descrittiva.
