Uma comparação bayesiana dos resultados do abandono escolar com R e brms

O conjunto de dados tem 14 campos.

• Código VCAA – ID administrativo para cada código
• Nome da escola
• Sector – G = Governo, C = Católica & I = Independente
• Localidade ou subúrbio
• Total de alunos que concluem o 12º ano
• Consentidores no caminho
• Requeridos
• O resto das colunas representa a percentagem de inquiridos para cada resultado

Para efeitos do nosso estudo transversal, estamos interessados nas proporções de resultados por sector, pelo que temos de converter este quadro de dados alargado num formato longo.

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)

Análise exploratória de dados

Vamos visualizar e resumir brevemente o nosso conjunto de dados. O sector governamental tem 157, o sector independente tem 57 e o sector católico tem 50 escolas listadas neste conjunto de dados.

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()

As distribuições contam uma história interessante. O resultado da licenciatura mostra a maior variabilidade em todos os sectores, com as escolas independentes a terem a maior proporção média de alunos que optam por prosseguir os estudos universitários. É interessante notar que as escolas públicas têm a maior proporção média de alunos que vão para o emprego depois de concluírem o ensino secundário. Todos os outros resultados não variam muito – voltaremos a este assunto em breve.

Modelação estatística – Regressão de verosimilhança beta

Estamos focados nas proporções de alunos por escola e nos seus resultados após a conclusão do ensino secundário. A probabilidade beta é a nossa melhor aposta nestes casos. brms é um pacote brilhante da Buerkner para desenvolver modelos Bayesianos. O nosso objetivo é modelar a distribuição das proporções por resultado e por sector.

A regressão beta modela dois parâmetros, μ e φ. μ representa a proporção média e φ é uma medida de precisão (ou variância).

O nosso primeiro modelo está representado abaixo, note que já estamos a começar com uma interação entre Sector e Resultado porque esta é a pergunta a que queremos que o modelo responda e, por isso, este é um modelo ANOVA.

Estamos a pedir ao modelo que crie termos Beta individuais para cada combinação de Sector e Resultado, com um termo φ agrupado, ou que modele diferentes médias de proporção com a mesma variância. Esperamos que 50% destas proporções se situem entre (-3, 1) na escala logit ou (0,01, 0,73) como proporções. Trata-se de uma priorização suficientemente ampla, mas informada. A estimativa global de Phi é um número positivo, pelo que utilizamos uma distribuição gama que é suficientemente ampla.

# Modelling Proportion with Sector:Outcome Interaction term using Beta - Pooled Phi term
m3a <- 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)

Observe que, na configuração do modelo, adicionamos um incremento de 1% a todas as proporções – isso ocorre porque a regressão beta não lida com valores zero. Tentámos modelar isto com beta inflacionado zero, mas demorou muito mais tempo a convergir.

Da mesma forma, podemos modelar sem um phi agrupado, o que intuitivamente faz sentido dado o que observámos nas distribuições acima, há variação entre cada uma das combinações de sector e resultado. O modelo é definido a seguir.

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)

Diagnóstico e comparação de modelos

Com dois modelos em mãos, vamos agora comparar a sua precisão preditiva fora da amostra usando a estimativa Bayesiana LOO da densidade preditiva log pointwise esperada (elpd_loo).

comp <- loo_compare(m3a, m3b)
print(comp, simplify = F)

Em termos simples, quanto maior for o valor esperado do logpoint wise leave one out, maior será a precisão da previsão em dados não vistos. Isto dá-nos uma boa relativa medida relativa da exatidão do modelo entre modelos. Pode ainda verificar isto através de uma verificação preditiva posterior, uma comparação entre valores observados e simulados. No nosso caso, o modelo m3b é o que melhor modela os dados observados.

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')

O modelo m3b, dada a sua variância não agrupada ou termo phi, é capaz de captar melhor a variação nas distribuições de proporções por sector e resultado.

ANOVA – Estilo Bayesiano

Lembre-se de que a nossa questão de investigação consiste em tentar compreender se existem diferenças nos resultados entre sectores e em que medida. Na estatística frequentista, poderíamos utilizar a ANOVA, uma abordagem de diferença de médias entre grupos. O ponto fraco desta abordagem é o facto de os resultados fornecerem uma estimativa e um intervalo de confiança, sem que se perceba a incerteza destas estimativas, e um valor p contra-intuitivo que indica se a diferença de médias é ou não estatisticamente significativa. Não, obrigado.

Abaixo, geramos um conjunto de valores esperados para cada interação de combinação de sector e resultado e, em seguida, utilizamos a excelente função tidybayes::compare_levels() que faz o trabalho pesado. Calcula a diferença nas médias posteriores entre sectores para cada resultado. Excluímos o resultado “outro” por uma questão de brevidade.

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')

Em alternativa, podemos resumir todas estas distribuições numa tabela simples para uma interpretação mais fácil e um intervalo de credibilidade de 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()

Por exemplo, se tivéssemos de resumir uma comparação num relatório formal, poderíamos escrever o seguinte.

Os alunos das escolas públicas têm menos probabilidades de iniciar estudos universitários do que os seus homólogos das escolas católicas e independentes após a conclusão do VCE.

Em média, 42,5% (entre 39,5% e 45,6%) dos alunos das escolas públicas, 53,2% (entre 47,7% e 58,4%) dos alunos das escolas católicas e 65% (entre 60,1% e 69,7%) dos alunos das escolas independentes iniciam o ensino superior após a conclusão do 12.

Existe uma probabilidade a posteriori de 89% de que a diferença entre as matrículas de alunos de escolas católicas e públicas se situe entre 5,6% e 15,7%, com uma média de 10,7%. Além disso, a diferença entre as matrículas de estudantes independentes e do Governo situa-se entre 17,8% e 27%, com uma média de 22,5%.

Estas diferenças são substanciais e existe uma probabilidade de 100% de que não sejam nulas.

Resumo e conclusão

Neste artigo demonstrámos como modelar dados proporcionais usando uma função de verosimilhança Beta usando modelação Bayesiana, e depois executar ANOVA Bayesiana para explorar as diferenças nos resultados proporcionais entre sectores.

Não procurámos criar uma compreensão causal destas diferenças. É possível imaginar que existem vários factores que influenciam o desempenho de cada aluno, o contexto socioeconómico, o nível de escolaridade dos pais, para além dos impactos a nível da escola, dos recursos, etc. Trata-se de uma área incrivelmente complexa da política pública que pode muitas vezes ficar atolada numa retórica de soma zero.

Pessoalmente, sou a primeira pessoa da minha família alargada a frequentar e a concluir o ensino superior. Frequentei uma escola secundária pública de nível médio e obtive resultados razoavelmente bons para entrar na minha primeira opção. Os meus pais encorajaram-me a prosseguir os estudos, tendo ambos optado por abandonar a escola aos 16 anos. Embora este artigo forneça provas de que a diferença entre escolas públicas e não públicas é real, tem um carácter puramente descritivo.